diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 3f240d9b52..9fa73efbfa 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -251,7 +251,6 @@ OBJS_ESOLVER=esolver.o\ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\ esolver_ks_lcao_tddft.o\ dpks_cal_e_delta_band.o\ - set_matrix_grid.o\ lcao_before_scf.o\ esolver_gets.o\ lcao_others.o\ diff --git a/source/module_esolver/CMakeLists.txt b/source/module_esolver/CMakeLists.txt index 77b8b0e9cc..1354c35707 100644 --- a/source/module_esolver/CMakeLists.txt +++ b/source/module_esolver/CMakeLists.txt @@ -18,7 +18,6 @@ if(ENABLE_LCAO) esolver_ks_lcao.cpp esolver_ks_lcao_tddft.cpp dpks_cal_e_delta_band.cpp - set_matrix_grid.cpp lcao_before_scf.cpp esolver_gets.cpp lcao_others.cpp diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index ea4bb2f879..0c376bc421 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -39,13 +39,94 @@ namespace ModuleESolver { template -void ESolver_KS_LCAO::beforesolver(const int istep) +void ESolver_KS_LCAO::before_scf(const int istep) { - ModuleBase::TITLE("ESolver_KS_LCAO", "beforesolver"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); + ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf"); + + //! 1) call before_scf() of ESolver_FP + ESolver_FP::before_scf(istep); + + if (GlobalC::ucell.ionic_position_updated) + { + this->CE.update_all_dis(GlobalC::ucell); + this->CE.extrapolate_charge( +#ifdef __MPI + &(GlobalC::Pgrid), +#endif + GlobalC::ucell, + this->pelec->charge, + &(this->sf), + GlobalV::ofs_running, + GlobalV::ofs_warning); + } + + //---------------------------------------------------------- + // about vdw, jiyy add vdwd3 and linpz add vdwd2 + //---------------------------------------------------------- + auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running)); + if (vdw_solver != nullptr) + { + this->pelec->f_en.evdw = vdw_solver->get_energy(); + } // 1. prepare HS matrices, prepare grid integral - this->set_matrix_grid(this->RA); + // (1) Find adjacent atoms for each atom. + double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running, + PARAM.inp.out_level, + orb_.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + PARAM.globalv.gamma_only_local); + + atom_arrange::search(PARAM.inp.search_pbc, + GlobalV::ofs_running, + GlobalC::GridD, + GlobalC::ucell, + search_radius, + PARAM.inp.test_atom_input); + + // (3) Periodic condition search for each grid. + double dr_uniform = 0.001; + std::vector rcuts; + std::vector> psi_u; + std::vector> dpsi_u; + std::vector> d2psi_u; + + Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb_, psi_u, dpsi_u, d2psi_u); + + this->GridT.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, + GlobalC::ucell, + GlobalC::GridD, + dr_uniform, + rcuts, + psi_u, + dpsi_u, + d2psi_u, + PARAM.inp.nstream); + psi_u.clear(); + psi_u.shrink_to_fit(); + dpsi_u.clear(); + dpsi_u.shrink_to_fit(); + d2psi_u.clear(); + d2psi_u.shrink_to_fit(); + + // (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(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation @@ -63,12 +144,10 @@ void ESolver_KS_LCAO::beforesolver(const int istep) { nsk = PARAM.inp.nspin; ncol = this->pv.ncol_bands; - if (PARAM.inp.ks_solver == "genelpa" - || PARAM.inp.ks_solver == "elpa" - || PARAM.inp.ks_solver == "lapack" - || PARAM.inp.ks_solver == "pexsi" - || PARAM.inp.ks_solver == "cusolver" - || PARAM.inp.ks_solver == "cusolvermp") { + if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack" + || PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver" + || PARAM.inp.ks_solver == "cusolvermp") + { ncol = this->pv.ncol; } } @@ -85,12 +164,11 @@ void ESolver_KS_LCAO::beforesolver(const int istep) } // init wfc from file - if(istep == 0 && PARAM.inp.init_wfc == "file") + if (istep == 0 && PARAM.inp.init_wfc == "file") { - if (! ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec)) + if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec)) { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::beforesolver", - "read wfc nao failed"); + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::beforesolver", "read wfc nao failed"); } } @@ -116,10 +194,11 @@ void ESolver_KS_LCAO::beforesolver(const int istep) orb_, DM #ifdef __EXX - , istep - , GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step - , GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr - , GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs + , + istep, + GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step, + GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr, + GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs #endif ); } @@ -169,41 +248,6 @@ void ESolver_KS_LCAO::beforesolver(const int istep) { GlobalC::ucell.cal_ux(); } - ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); -} - -template -void ESolver_KS_LCAO::before_scf(const int istep) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf"); - - //! 1) call before_scf() of ESolver_FP - ESolver_FP::before_scf(istep); - - if (GlobalC::ucell.ionic_position_updated) - { - this->CE.update_all_dis(GlobalC::ucell); - this->CE.extrapolate_charge( -#ifdef __MPI - &(GlobalC::Pgrid), -#endif - GlobalC::ucell, - this->pelec->charge, - &(this->sf), - GlobalV::ofs_running, - GlobalV::ofs_warning); - } - - //---------------------------------------------------------- - // about vdw, jiyy add vdwd3 and linpz add vdwd2 - //---------------------------------------------------------- - auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running)); - if (vdw_solver != nullptr) - { - this->pelec->f_en.evdw = vdw_solver->get_energy(); - } - - this->beforesolver(istep); // Peize Lin add 2016-12-03 #ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 07c3518473..2d70c02611 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -70,7 +70,186 @@ void ESolver_KS_LCAO::others(const int istep) return; } - this->beforesolver(istep); + // 1. prepare HS matrices, prepare grid integral + // (1) Find adjacent atoms for each atom. + double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running, + PARAM.inp.out_level, + orb_.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + PARAM.globalv.gamma_only_local); + + atom_arrange::search(PARAM.inp.search_pbc, + GlobalV::ofs_running, + GlobalC::GridD, + GlobalC::ucell, + search_radius, + PARAM.inp.test_atom_input); + + // (3) Periodic condition search for each grid. + double dr_uniform = 0.001; + std::vector rcuts; + std::vector> psi_u; + std::vector> dpsi_u; + std::vector> d2psi_u; + + Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb_, psi_u, dpsi_u, d2psi_u); + + this->GridT.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, + GlobalC::ucell, + GlobalC::GridD, + dr_uniform, + rcuts, + psi_u, + dpsi_u, + d2psi_u, + PARAM.inp.nstream); + psi_u.clear(); + psi_u.shrink_to_fit(); + dpsi_u.clear(); + dpsi_u.shrink_to_fit(); + d2psi_u.clear(); + d2psi_u.shrink_to_fit(); + + // (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(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + + // 2. density matrix extrapolation + + // set the augmented orbitals index. + // after ParaO and GridT, + // this information is used to calculate + // the force. + + // init psi + if (this->psi == nullptr) + { + int nsk = 0; + int ncol = 0; + if (PARAM.globalv.gamma_only_local) + { + nsk = PARAM.inp.nspin; + ncol = this->pv.ncol_bands; + if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack" + || PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver" + || PARAM.inp.ks_solver == "cusolvermp") + { + ncol = this->pv.ncol; + } + } + else + { + nsk = this->kv.get_nks(); +#ifdef __MPI + ncol = this->pv.ncol_bands; +#else + ncol = PARAM.inp.nbands; +#endif + } + this->psi = new psi::Psi(nsk, ncol, this->pv.nrow, nullptr); + } + + // init wfc from file + if (istep == 0 && PARAM.inp.init_wfc == "file") + { + if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, this->pv, *(this->psi), this->pelec)) + { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::beforesolver", "read wfc nao failed"); + } + } + + // prepare grid in Gint + LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, orb_, *this->pw_rho, *this->pw_big); + + // init Hamiltonian + if (this->p_hamilt != nullptr) + { + delete this->p_hamilt; + this->p_hamilt = nullptr; + } + if (this->p_hamilt == nullptr) + { + elecstate::DensityMatrix* DM = dynamic_cast*>(this->pelec)->get_DM(); + this->p_hamilt = new hamilt::HamiltLCAO( + PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, + PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), + &this->pv, + this->pelec->pot, + this->kv, + two_center_bundle_, + orb_, + DM +#ifdef __EXX + , + istep, + GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step, + GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr, + GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs +#endif + ); + } + +#ifdef __DEEPKS + // for each ionic step, the overlap must be rebuilt + // since it depends on ionic positions + if (PARAM.globalv.deepks_setorb) + { + const Parallel_Orbitals* pv = &this->pv; + // build and save at beginning + GlobalC::ld.build_psialpha(PARAM.inp.cal_force, + GlobalC::ucell, + orb_, + GlobalC::GridD, + *(two_center_bundle_.overlap_orb_alpha)); + + if (PARAM.inp.deepks_out_unittest) + { + GlobalC::ld.check_psialpha(PARAM.inp.cal_force, GlobalC::ucell, orb_, GlobalC::GridD); + } + } +#endif + if (PARAM.inp.sc_mag_switch) + { + spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); + sc.init_sc(PARAM.inp.sc_thr, + PARAM.inp.nsc, + PARAM.inp.nsc_min, + PARAM.inp.alpha_trial, + PARAM.inp.sccut, + PARAM.inp.sc_drop_thr, + GlobalC::ucell, + &(this->pv), + PARAM.inp.nspin, + this->kv, + PARAM.inp.ks_solver, + this->p_hamilt, + this->psi, + this->pelec); + } + //========================================================= + // cal_ux should be called before init_scf because + // the direction of ux is used in noncoline_rho + //========================================================= + if (PARAM.inp.nspin == 4) + { + GlobalC::ucell.cal_ux(); + } + // pelec should be initialized before these calculations this->pelec->init_scf(istep, this->sf.strucFac, GlobalC::ucell.symm); // self consistent calculations for electronic ground state diff --git a/source/module_esolver/set_matrix_grid.cpp b/source/module_esolver/set_matrix_grid.cpp deleted file mode 100644 index 4a20d789eb..0000000000 --- a/source/module_esolver/set_matrix_grid.cpp +++ /dev/null @@ -1,98 +0,0 @@ -#include "esolver_ks_lcao.h" -#include "module_base/timer.h" -#include "module_cell/module_neighbor/sltk_atom_arrange.h" -#include "module_cell/module_neighbor/sltk_grid_driver.h" -#include "module_elecstate/elecstate_lcao.h" -#include "module_elecstate/module_charge/symmetry_rho.h" -#include "module_elecstate/module_dm/cal_dm_psi.h" -#include "module_esolver/esolver_ks_lcao.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" -#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" -#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" -#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_io/read_wfc_nao.h" -#include "module_io/write_elecstat_pot.h" -#include "module_io/write_wfc_nao.h" -#include "module_parameter/parameter.h" - -namespace ModuleESolver -{ - -template -void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "set_matrix_grid"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "set_matrix_grid"); - - // (1) Find adjacent atoms for each atom. - double search_radius = -1.0; - search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running, - PARAM.inp.out_level, - orb_.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), - PARAM.globalv.gamma_only_local); - - atom_arrange::search(PARAM.inp.search_pbc, - GlobalV::ofs_running, - GlobalC::GridD, - GlobalC::ucell, - search_radius, - PARAM.inp.test_atom_input); - - // ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"SEARCH ADJACENT - // ATOMS"); - - // (3) Periodic condition search for each grid. - double dr_uniform = 0.001; - std::vector rcuts; - std::vector> psi_u; - std::vector> dpsi_u; - std::vector> d2psi_u; - - Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb_, psi_u, dpsi_u, d2psi_u); - - this->GridT.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, - GlobalC::ucell, - GlobalC::GridD, - dr_uniform, - rcuts, - psi_u, - dpsi_u, - d2psi_u, - PARAM.inp.nstream); - psi_u.clear(); - psi_u.shrink_to_fit(); - dpsi_u.clear(); - dpsi_u.shrink_to_fit(); - d2psi_u.clear(); - d2psi_u.shrink_to_fit(); - - // (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. - ra.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); - - ModuleBase::timer::tick("ESolver_KS_LCAO", "set_matrix_grid"); - return; -} - -template class ESolver_KS_LCAO; -template class ESolver_KS_LCAO, double>; -template class ESolver_KS_LCAO, std::complex>; - -}