diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 19948c6a222..537b631ffa4 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -266,6 +266,7 @@ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\ lcao_after_scf.o\ esolver_gets.o\ lcao_others.o\ + esolver_dm2rho.o\ OBJS_GINT=gint.o\ gint_gamma_env.o\ diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index 024fa3e45c3..efa7b1d76bd 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -23,32 +23,6 @@ void ElecStateLCAO>::psiToRho(const psi::Psi rho are, technically, two separate - // // functionalities, as we cannot rule out the possibility that we may have a - // // dm from other sources, such as read from file. However, since we are not - // // separating them now, I opt to add a flag to control how dm is obtained as - // // of now - // if (!PARAM.inp.dm_to_rho) - // { - // ModuleBase::GlobalFunc::NOTE("Calculate the density matrix."); - - // // this part for calculating DMK in 2d-block format, not used for charge - // // now - // // psi::Psi> dm_k_2d(); - - // if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == - // "scalapack_gvx" || PARAM.inp.ks_solver == "lapack" - // || PARAM.inp.ks_solver == "cusolver" || PARAM.inp.ks_solver == "cusolvermp" - // || PARAM.inp.ks_solver == "cg_in_lcao") // Peize Lin test 2019-05-15 - // { - // elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), - // this->wg, - // psi, - // *(this->DM)); - // this->DM->cal_DMR(); - // } - // } - for (int is = 0; is < PARAM.inp.nspin; is++) { ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], diff --git a/source/module_esolver/CMakeLists.txt b/source/module_esolver/CMakeLists.txt index 6b76e8f84ed..c9fda449482 100644 --- a/source/module_esolver/CMakeLists.txt +++ b/source/module_esolver/CMakeLists.txt @@ -20,6 +20,7 @@ if(ENABLE_LCAO) lcao_after_scf.cpp esolver_gets.cpp lcao_others.cpp + esolver_dm2rho.cpp ) endif() diff --git a/source/module_esolver/esolver.cpp b/source/module_esolver/esolver.cpp index 46f4e238401..9d5559dd0d7 100644 --- a/source/module_esolver/esolver.cpp +++ b/source/module_esolver/esolver.cpp @@ -5,6 +5,7 @@ #include "module_base/module_device/device.h" #include "module_parameter/parameter.h" #ifdef __LCAO +#include "esolver_dm2rho.h" #include "esolver_gets.h" #include "esolver_ks_lcao.h" #include "esolver_ks_lcao_tddft.h" @@ -204,11 +205,25 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell) } else if (PARAM.inp.nspin < 4) { - return new ESolver_KS_LCAO, double>(); + if (PARAM.inp.dm_to_rho) + { + return new ESolver_DM2rho, double>(); + } + else + { + return new ESolver_KS_LCAO, double>(); + } } else { - return new ESolver_KS_LCAO, std::complex>(); + if (PARAM.inp.dm_to_rho) + { + return new ESolver_DM2rho, std::complex>(); + } + else + { + return new ESolver_KS_LCAO, std::complex>(); + } } } else if (esolver_type == "ksdft_lcao_tddft") diff --git a/source/module_esolver/esolver_dm2rho.cpp b/source/module_esolver/esolver_dm2rho.cpp new file mode 100644 index 00000000000..2a10c713067 --- /dev/null +++ b/source/module_esolver/esolver_dm2rho.cpp @@ -0,0 +1,100 @@ +#include "esolver_dm2rho.h" + +#include "module_base/timer.h" +#include "module_cell/module_neighbor/sltk_atom_arrange.h" +#include "module_elecstate/elecstate_lcao.h" +#include "module_elecstate/read_pseudo.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/operator_lcao.h" +#include "module_io/cube_io.h" +#include "module_io/io_npz.h" +#include "module_io/print_info.h" + +namespace ModuleESolver +{ + +template +ESolver_DM2rho::ESolver_DM2rho() +{ + this->classname = "ESolver_DM2rho"; + this->basisname = "LCAO"; +} + +template +ESolver_DM2rho::~ESolver_DM2rho() +{ +} + +template +void ESolver_DM2rho::before_all_runners(UnitCell& ucell, const Input_para& inp) +{ + ModuleBase::TITLE("ESolver_DM2rho", "before_all_runners"); + ModuleBase::timer::tick("ESolver_DM2rho", "before_all_runners"); + + ESolver_KS_LCAO::before_all_runners(ucell, inp); + + ModuleBase::timer::tick("ESolver_DM2rho", "before_all_runners"); +} + +template +void ESolver_DM2rho::runner(UnitCell& ucell, const int istep) +{ + ModuleBase::TITLE("ESolver_DM2rho", "runner"); + ModuleBase::timer::tick("ESolver_DM2rho", "runner"); + + ESolver_KS_LCAO::before_scf(ucell, istep); + + // file name of DM + std::string zipname = "output_DM0.npz"; + elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); + + // read DM from file + ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1))); + + // if nspin=2, need extra reading + if (PARAM.inp.nspin == 2) + { + zipname = "output_DM1.npz"; + ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2))); + } + + this->pelec->psiToRho(*this->psi); + + int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1; + + for (int is = 0; is < nspin0; is++) + { + std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube"; + + // write electron density + ModuleIO::write_vdata_palgrid(this->Pgrid, + this->chr.rho[is], + is, + PARAM.inp.nspin, + istep, + fn, + this->pelec->eferm.get_efval(is), + &(ucell), + 3, + 1); + } + + ModuleBase::timer::tick("ESolver_DM2rho", "runner"); +} + +template +void ESolver_DM2rho::after_all_runners(UnitCell& ucell) +{ + ModuleBase::TITLE("ESolver_DM2rho", "after_all_runners"); + ModuleBase::timer::tick("ESolver_DM2rho", "after_all_runners"); + + ESolver_KS_LCAO::after_all_runners(ucell); + + ModuleBase::timer::tick("ESolver_DM2rho", "after_all_runners"); +}; + +template class ESolver_DM2rho, double>; +template class ESolver_DM2rho, std::complex>; + +} // namespace ModuleESolver diff --git a/source/module_esolver/esolver_dm2rho.h b/source/module_esolver/esolver_dm2rho.h new file mode 100644 index 00000000000..1e28ac127fb --- /dev/null +++ b/source/module_esolver/esolver_dm2rho.h @@ -0,0 +1,25 @@ +#ifndef ESOLVER_DM2RHO_H +#define ESOLVER_DM2RHO_H + +#include "module_esolver/esolver_ks_lcao.h" + +#include + +namespace ModuleESolver +{ + +template +class ESolver_DM2rho : public ESolver_KS_LCAO +{ + public: + ESolver_DM2rho(); + ~ESolver_DM2rho(); + + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; + + void after_all_runners(UnitCell& ucell) override; + + void runner(UnitCell& ucell, const int istep) override; +}; +} // namespace ModuleESolver +#endif diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 1f31453ebda..f99fb2264b7 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -152,20 +152,10 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso { for (int is = 0; is < PARAM.inp.nspin; is++) { - double* data = nullptr; - if (PARAM.inp.dm_to_rho) - { - data = this->chr.rho[is]; - this->pw_rhod->real2recip(this->chr.rho[is], this->chr.rhog[is]); - } - else - { - data = this->chr.rho_save[is]; - this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]); - } + this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]); std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube"; ModuleIO::write_vdata_palgrid(Pgrid, - data, + this->chr.rho_save[is], is, PARAM.inp.nspin, istep, @@ -381,19 +371,16 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& { if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver) { - std::complex** rhog_tot - = (PARAM.inp.dm_to_rho) ? this->chr.rhog : this->chr.rhog_save; - double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->chr.rho : this->chr.rho_save; for (int is = 0; is < PARAM.inp.nspin; is++) { - this->pw_rhod->real2recip(rhor_tot[is], rhog_tot[is]); + this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]); } ModuleIO::write_rhog(PARAM.globalv.global_out_dir + PARAM.inp.suffix + "-CHARGE-DENSITY.restart", PARAM.globalv.gamma_only_pw || PARAM.globalv.gamma_only_local, this->pw_rhod, PARAM.inp.nspin, ucell.GT, - rhog_tot, + this->chr.rhog_save, GlobalV::MY_POOL, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL); diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 932da275293..2a26b2c7de7 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -439,13 +439,6 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) // 2) before_scf (electronic iteration loops) this->before_scf(ucell, istep); - // 3) write charge density - if (PARAM.inp.dm_to_rho) - { - ModuleBase::timer::tick(this->classname, "runner"); - return; // nothing further is needed - } - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF"); // 4) SCF iterations diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 9549cc511e9..7020e2e9fb5 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -637,13 +637,11 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const this->pelec->nelec_spin, this->pelec->skip_weights); - if (!PARAM.inp.dm_to_rho) - { - auto _pelec = dynamic_cast*>(this->pelec); - elecstate::calEBand(_pelec->ekb,_pelec->wg,_pelec->f_en); - elecstate::cal_dm_psi(_pelec->DM->get_paraV_pointer(), _pelec->wg, *this->psi, *(_pelec->DM)); - _pelec->DM->cal_DMR(); - } + auto _pelec = dynamic_cast*>(this->pelec); + elecstate::calEBand(_pelec->ekb, _pelec->wg, _pelec->f_en); + elecstate::cal_dm_psi(_pelec->DM->get_paraV_pointer(), _pelec->wg, *this->psi, *(_pelec->DM)); + _pelec->DM->cal_DMR(); + this->pelec->psiToRho(*this->psi); this->pelec->skip_weights = false; diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 0051a50c191..08b3e47c556 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -251,57 +251,6 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) dynamic_cast*>(this->pelec)->get_DM()->cal_DMR(); } - if (PARAM.inp.dm_to_rho) - { - // file name of DM - std::string zipname = "output_DM0.npz"; - elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); - - // read DM from file - ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1))); - - // if nspin=2, need extra reading - if (PARAM.inp.nspin == 2) - { - zipname = "output_DM1.npz"; - ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2))); - } - - elecstate::calculate_weights(this->pelec->ekb, - this->pelec->wg, - this->pelec->klist, - this->pelec->eferm, - this->pelec->f_en, - this->pelec->nelec_spin, - this->pelec->skip_weights); - - this->pelec->psiToRho(*this->psi); - - int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1; - - for (int is = 0; is < nspin0; is++) - { - std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube"; - - // write electron density - ModuleIO::write_vdata_palgrid(this->Pgrid, - this->chr.rho[is], - is, - PARAM.inp.nspin, - istep, - fn, - this->pelec->eferm.get_efval(is), - &(ucell), - 3, - 1); - } - - // why we need to return here? mohan add 2025-03-10 - ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); - return; - } - // 16) the electron charge density should be symmetrized, // here is the initialization Symmetry_rho srho; diff --git a/source/module_hsolver/hsolver_lcao.cpp b/source/module_hsolver/hsolver_lcao.cpp index 1e0cfa6accc..6a361e08db6 100644 --- a/source/module_hsolver/hsolver_lcao.cpp +++ b/source/module_hsolver/hsolver_lcao.cpp @@ -83,13 +83,11 @@ void HSolverLCAO::solve(hamilt::Hamilt* pHamilt, pes->f_en, pes->nelec_spin, pes->skip_weights); - if (!PARAM.inp.dm_to_rho) - { - auto _pes_lcao = dynamic_cast*>(pes); - elecstate::calEBand(_pes_lcao->ekb,_pes_lcao->wg,_pes_lcao->f_en); - elecstate::cal_dm_psi(_pes_lcao->DM->get_paraV_pointer(), _pes_lcao->wg, psi, *(_pes_lcao->DM)); - _pes_lcao->DM->cal_DMR(); - } + + auto _pes_lcao = dynamic_cast*>(pes); + elecstate::calEBand(_pes_lcao->ekb, _pes_lcao->wg, _pes_lcao->f_en); + elecstate::cal_dm_psi(_pes_lcao->DM->get_paraV_pointer(), _pes_lcao->wg, psi, *(_pes_lcao->DM)); + _pes_lcao->DM->cal_DMR(); if (!skip_charge) {