From c03c8796f382b89a45d29c550c06197a4472d137 Mon Sep 17 00:00:00 2001 From: Yu Liu <77716030+YuLiu98@users.noreply.github.com> Date: Mon, 11 Nov 2024 16:33:04 +0800 Subject: [PATCH] Refactor: runner() of esolver_ks (#5445) * Refactor: runner() of esolver_ks * rename hamilt2density and diag as hamilt2density_single and hamilt2density --- source/module_esolver/esolver_ks.cpp | 324 ++++++++++-------- source/module_esolver/esolver_ks.h | 8 +- source/module_esolver/esolver_ks_lcao.cpp | 137 +++----- source/module_esolver/esolver_ks_lcao.h | 4 +- .../module_esolver/esolver_ks_lcao_tddft.cpp | 51 +-- source/module_esolver/esolver_ks_lcao_tddft.h | 4 +- source/module_esolver/esolver_ks_lcaopw.cpp | 86 ++--- source/module_esolver/esolver_ks_lcaopw.h | 2 +- source/module_esolver/esolver_ks_pw.cpp | 111 ++---- source/module_esolver/esolver_ks_pw.h | 2 +- source/module_esolver/esolver_sdft_pw.cpp | 4 +- source/module_esolver/esolver_sdft_pw.h | 2 +- 12 files changed, 332 insertions(+), 403 deletions(-) diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 4d3121a904..8f9a715d5a 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -365,18 +365,72 @@ void ESolver_KS::init_after_vc(const Input_para& inp, UnitCell& ucell } //------------------------------------------------------------------------------ -//! the 5th function of ESolver_KS: hamilt2density +//! the 5th function of ESolver_KS: hamilt2density_single //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS::hamilt2density(const int istep, const int iter, const double ethr) +void ESolver_KS::hamilt2density_single(const int istep, const int iter, const double ethr) { - ModuleBase::timer::tick(this->classname, "hamilt2density"); + ModuleBase::timer::tick(this->classname, "hamilt2density_single"); // Temporarily, before HSolver is constructed, it should be overrided by // LCAO, PW, SDFT and TDDFT. // After HSolver is constructed, LCAO, PW, SDFT should delete their own - // hamilt2density() and use: - ModuleBase::timer::tick(this->classname, "hamilt2density"); + // hamilt2density_single() and use: + ModuleBase::timer::tick(this->classname, "hamilt2density_single"); +} + +template +void ESolver_KS::hamilt2density(const int istep, const int iter, const double ethr) +{ + // 7) use Hamiltonian to obtain charge density + this->hamilt2density_single(istep, iter, diag_ethr); + + // 8) for MPI: STOGROUP? need to rewrite + // It may be changed when more clever parallel algorithm is + // put forward. + // When parallel algorithm for bands are adopted. Density will only be + // treated in the first group. + //(Different ranks should have abtained the same, but small differences + // always exist in practice.) + // Maybe in the future, density and wavefunctions should use different + // parallel algorithms, in which they do not occupy all processors, for + // example wavefunctions uses 20 processors while density uses 10. + if (GlobalV::MY_STOGROUP == 0) + { + // double drho = this->estate.caldr2(); + // EState should be used after it is constructed. + + drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); + hsolver_error = 0.0; + if (iter == 1) + { + hsolver_error + = hsolver::cal_hsolve_error(PARAM.inp.basis_type, PARAM.inp.esolver_type, diag_ethr, PARAM.inp.nelec); + + // The error of HSolver is larger than drho, + // so a more precise HSolver should be executed. + if (hsolver_error > drho) + { + diag_ethr = hsolver::reset_diag_ethr(GlobalV::ofs_running, + PARAM.inp.basis_type, + PARAM.inp.esolver_type, + PARAM.inp.precision, + hsolver_error, + drho, + diag_ethr, + PARAM.inp.nelec); + + this->hamilt2density_single(istep, iter, diag_ethr); + + drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); + + hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type, + PARAM.inp.esolver_type, + diag_ethr, + PARAM.inp.nelec); + } + } + } } //------------------------------------------------------------------------------ @@ -418,7 +472,6 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF"); - bool firstscf = true; this->conv_esolver = false; this->niter = this->maxniter; @@ -431,146 +484,8 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) // 6) initialization of SCF iterations this->iter_init(istep, iter); - // 7) use Hamiltonian to obtain charge density this->hamilt2density(istep, iter, diag_ethr); - // 8) for MPI: STOGROUP? need to rewrite - // It may be changed when more clever parallel algorithm is - // put forward. - // When parallel algorithm for bands are adopted. Density will only be - // treated in the first group. - //(Different ranks should have abtained the same, but small differences - // always exist in practice.) - // Maybe in the future, density and wavefunctions should use different - // parallel algorithms, in which they do not occupy all processors, for - // example wavefunctions uses 20 processors while density uses 10. - if (GlobalV::MY_STOGROUP == 0) - { - // double drho = this->estate.caldr2(); - // EState should be used after it is constructed. - - drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); - double hsolver_error = 0.0; - if (firstscf) - { - firstscf = false; - hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type, - PARAM.inp.esolver_type, - diag_ethr, - PARAM.inp.nelec); - - // The error of HSolver is larger than drho, - // so a more precise HSolver should be executed. - if (hsolver_error > drho) - { - diag_ethr = hsolver::reset_diag_ethr(GlobalV::ofs_running, - PARAM.inp.basis_type, - PARAM.inp.esolver_type, - PARAM.inp.precision, - hsolver_error, - drho, - diag_ethr, - PARAM.inp.nelec); - - this->hamilt2density(istep, iter, diag_ethr); - - drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); - - hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type, - PARAM.inp.esolver_type, - diag_ethr, - PARAM.inp.nelec); - } - } - // mixing will restart at this->p_chgmix->mixing_restart steps - if (drho <= PARAM.inp.mixing_restart && PARAM.inp.mixing_restart > 0.0 - && this->p_chgmix->mixing_restart_step > iter) - { - this->p_chgmix->mixing_restart_step = iter + 1; - } - - if (PARAM.inp.scf_os_stop) // if oscillation is detected, SCF will stop - { - this->oscillate_esolver = this->p_chgmix->if_scf_oscillate(iter, drho, PARAM.inp.scf_os_ndim, PARAM.inp.scf_os_thr); - } - - // drho will be 0 at this->p_chgmix->mixing_restart step, which is - // not ground state - bool not_restart_step = !(iter == this->p_chgmix->mixing_restart_step && PARAM.inp.mixing_restart > 0.0); - // SCF will continue if U is not converged for uramping calculation - bool is_U_converged = true; - // to avoid unnecessary dependence on dft+u, refactor is needed -#ifdef __LCAO - if (PARAM.inp.dft_plus_u) - { - is_U_converged = GlobalC::dftu.u_converged(); - } -#endif - - this->conv_esolver = (drho < this->scf_thr && not_restart_step && is_U_converged); - - // add energy threshold for SCF convergence - if (this->scf_ene_thr > 0.0) - { - // calculate energy of output charge density - this->update_pot(istep, iter); - this->pelec->cal_energies(2); // 2 means Kohn-Sham functional - // now, etot_old is the energy of input density, while etot is the energy of output density - this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old; - // output etot_delta - GlobalV::ofs_running << " DeltaE_womix = " << this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV << " eV" << std::endl; - if (iter > 1 && this->conv_esolver == 1) // only check when density is converged - { - // update the convergence flag - this->conv_esolver - = (std::abs(this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < this->scf_ene_thr); - } - } - - // If drho < hsolver_error in the first iter or drho < scf_thr, we - // do not change rho. - if (drho < hsolver_error || this->conv_esolver) - { - if (drho < hsolver_error) - { - GlobalV::ofs_warning << " drho < hsolver_error, keep " - "charge density unchanged." - << std::endl; - } - } - else - { - //----------charge mixing--------------- - // mixing will restart after this->p_chgmix->mixing_restart - // steps - if (PARAM.inp.mixing_restart > 0 && iter == this->p_chgmix->mixing_restart_step - 1 - && drho <= PARAM.inp.mixing_restart) - { - // do not mix charge density - } - else - { - p_chgmix->mix_rho(pelec->charge); // update chr->rho by mixing - } - if (PARAM.inp.scf_thr_type == 2) - { - pelec->charge->renormalize_rho(); // renormalize rho in R-space would - // induce a error in K-space - } - //----------charge mixing done----------- - } - } -#ifdef __MPI - MPI_Bcast(&drho, 1, MPI_DOUBLE, 0, PARAPW_WORLD); - MPI_Bcast(&this->conv_esolver, 1, MPI_DOUBLE, 0, PARAPW_WORLD); - MPI_Bcast(pelec->charge->rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, PARAPW_WORLD); -#endif - - // 9) update potential - // Hamilt should be used after it is constructed. - // this->phamilt->update(conv_esolver); - this->update_pot(istep, iter); - // 10) finish scf iterations this->iter_finish(istep, iter); @@ -635,13 +550,134 @@ void ESolver_KS::iter_init(const int istep, const int iter) PARAM.inp.nbands, esolver_KS_ne); } + + // 1) save input rho + this->pelec->charge->save_rho_before_sum_band(); } template void ESolver_KS::iter_finish(const int istep, int& iter) { + if (PARAM.inp.out_bandgap) + { + if (!PARAM.globalv.two_fermi) + { + this->pelec->cal_bandgap(); + } + else + { + this->pelec->cal_bandgap_updw(); + } + } + + for (int ik = 0; ik < this->kv.get_nks(); ++ik) + { + this->pelec->print_band(ik, PARAM.inp.printe, iter); + } + + // compute magnetization, only for LSDA(spin==2) + GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, + this->pelec->charge->nxyz, + this->pelec->charge->rho, + this->pelec->nelec_spin.data()); + + if (GlobalV::MY_STOGROUP == 0) + { + // mixing will restart at this->p_chgmix->mixing_restart steps + if (drho <= PARAM.inp.mixing_restart && PARAM.inp.mixing_restart > 0.0 + && this->p_chgmix->mixing_restart_step > iter) + { + this->p_chgmix->mixing_restart_step = iter + 1; + } + + if (PARAM.inp.scf_os_stop) // if oscillation is detected, SCF will stop + { + this->oscillate_esolver + = this->p_chgmix->if_scf_oscillate(iter, drho, PARAM.inp.scf_os_ndim, PARAM.inp.scf_os_thr); + } + + // drho will be 0 at this->p_chgmix->mixing_restart step, which is + // not ground state + bool not_restart_step = !(iter == this->p_chgmix->mixing_restart_step && PARAM.inp.mixing_restart > 0.0); + // SCF will continue if U is not converged for uramping calculation + bool is_U_converged = true; + // to avoid unnecessary dependence on dft+u, refactor is needed +#ifdef __LCAO + if (PARAM.inp.dft_plus_u) + { + is_U_converged = GlobalC::dftu.u_converged(); + } +#endif + + this->conv_esolver = (drho < this->scf_thr && not_restart_step && is_U_converged); + + // add energy threshold for SCF convergence + if (this->scf_ene_thr > 0.0) + { + // calculate energy of output charge density + this->update_pot(istep, iter); + this->pelec->cal_energies(2); // 2 means Kohn-Sham functional + // now, etot_old is the energy of input density, while etot is the energy of output density + this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old; + // output etot_delta + GlobalV::ofs_running << " DeltaE_womix = " << this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV << " eV" + << std::endl; + if (iter > 1 && this->conv_esolver == 1) // only check when density is converged + { + // update the convergence flag + this->conv_esolver + = (std::abs(this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < this->scf_ene_thr); + } + } + + // If drho < hsolver_error in the first iter or drho < scf_thr, we + // do not change rho. + if (drho < hsolver_error || this->conv_esolver) + { + if (drho < hsolver_error) + { + GlobalV::ofs_warning << " drho < hsolver_error, keep " + "charge density unchanged." + << std::endl; + } + } + else + { + //----------charge mixing--------------- + // mixing will restart after this->p_chgmix->mixing_restart + // steps + if (PARAM.inp.mixing_restart > 0 && iter == this->p_chgmix->mixing_restart_step - 1 + && drho <= PARAM.inp.mixing_restart) + { + // do not mix charge density + } + else + { + p_chgmix->mix_rho(pelec->charge); // update chr->rho by mixing + } + if (PARAM.inp.scf_thr_type == 2) + { + pelec->charge->renormalize_rho(); // renormalize rho in R-space would + // induce a error in K-space + } + //----------charge mixing done----------- + } + } + +#ifdef __MPI + MPI_Bcast(&drho, 1, MPI_DOUBLE, 0, PARAPW_WORLD); + MPI_Bcast(&this->conv_esolver, 1, MPI_DOUBLE, 0, PARAPW_WORLD); + MPI_Bcast(pelec->charge->rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, PARAPW_WORLD); +#endif + + // update potential + // Hamilt should be used after it is constructed. + // this->phamilt->update(conv_esolver); + this->update_pot(istep, iter); + // 1 means Harris-Foulkes functional // 2 means Kohn-Sham functional + this->pelec->cal_energies(1); this->pelec->cal_energies(2); if (iter == 1) diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 8b952a6e22..88d127b619 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -46,12 +46,15 @@ class ESolver_KS : public ESolver_FP //! Something to do after hamilt2density function in each iter loop. virtual void iter_finish(const int istep, int& iter); - // calculate electron density from a specific Hamiltonian - virtual void hamilt2density(const int istep, const int iter, const double ethr); + // calculate electron density from a specific Hamiltonian with ethr + virtual void hamilt2density_single(const int istep, const int iter, const double ethr); // calculate electron states from a specific Hamiltonian virtual void hamilt2estates(const double ethr) {}; + // calculate electron density from a specific Hamiltonian + void hamilt2density(const int istep, const int iter, const double ethr); + //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. virtual void after_scf(const int istep) override; @@ -82,6 +85,7 @@ class ESolver_KS : public ESolver_FP double scf_thr; // scf density threshold double scf_ene_thr; // scf energy threshold double drho; // the difference between rho_in (before HSolver) and rho_out (After HSolver) + double hsolver_error; // the error of HSolver int maxniter; // maximum iter steps for scf int niter; // iter steps actually used in scf int out_freq_elec; // frequency for output diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 5bd2dca2f4..0249f7e3ae 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -681,10 +681,17 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) SpinConstrain& sc = SpinConstrain::getScInstance(); sc.run_lambda_loop(iter - 1); } + + // save density matrix DMR for mixing + if (PARAM.inp.mixing_restart > 0 && PARAM.inp.mixing_dmr && this->p_chgmix->mixing_restart_count > 0) + { + elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); + dm->save_DMR(); + } } //------------------------------------------------------------------------------ -//! the 11th function of ESolver_KS_LCAO: hamilt2density +//! the 11th function of ESolver_KS_LCAO: hamilt2density_single //! mohan add 2024-05-11 //! 1) save input rho //! 2) save density matrix DMR for mixing @@ -700,47 +707,16 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) //! 12) calculate delta energy //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) +void ESolver_KS_LCAO::hamilt2density_single(int istep, int iter, double ethr) { - ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density"); - - // 1) save input rho - this->pelec->charge->save_rho_before_sum_band(); - - // 2) save density matrix DMR for mixing - if (PARAM.inp.mixing_restart > 0 && PARAM.inp.mixing_dmr && this->p_chgmix->mixing_restart_count > 0) - { - elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); - dm->save_DMR(); - } + ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density_single"); - // 3) solve the Hamiltonian and output band gap - { - // reset energy - this->pelec->f_en.eband = 0.0; - this->pelec->f_en.demet = 0.0; - - hsolver::HSolverLCAO hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver); - hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, false); + // reset energy + this->pelec->f_en.eband = 0.0; + this->pelec->f_en.demet = 0.0; - if (PARAM.inp.out_bandgap) - { - if (!PARAM.globalv.two_fermi) - { - this->pelec->cal_bandgap(); - } - else - { - this->pelec->cal_bandgap_updw(); - } - } - } - - // 4) print bands for each k-point and each band - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { - this->pelec->print_band(ik, PARAM.inp.printe, iter); - } + hsolver::HSolverLCAO hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver); + hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, false); // 5) what's the exd used for? #ifdef __EXX @@ -754,46 +730,6 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) } #endif - // 6) calculate the local occupation number matrix and energy correction in - // DFT+U - if (PARAM.inp.dft_plus_u) - { - // only old DFT+U method should calculated energy correction in esolver, - // new DFT+U method will calculate energy in calculating Hamiltonian - if (PARAM.inp.dft_plus_u == 2) - { - if (GlobalC::dftu.omc != 2) - { - const std::vector>& tmp_dm - = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); - this->dftu_cal_occup_m(iter, tmp_dm); - } - GlobalC::dftu.cal_energy_correction(istep); - } - GlobalC::dftu.output(); - } - - // (7) for deepks, calculate delta_e -#ifdef __DEEPKS - if (PARAM.inp.deepks_scf) - { - const std::vector>& dm - = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); - - this->dpks_cal_e_delta_band(dm); - } -#endif - - // 8) for delta spin - if (PARAM.inp.sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); - sc.cal_MW(iter, this->p_hamilt); - } - - // 9) use new charge density to calculate energy - this->pelec->cal_energies(1); - // 10) symmetrize the charge density Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) @@ -801,12 +737,6 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::ucell.symm); } - // 11) compute magnetization, only for spin==2 - GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, - this->pelec->charge->nxyz, - this->pelec->charge->rho, - this->pelec->nelec_spin.data()); - // 12) calculate delta energy this->pelec->f_en.deband = this->pelec->cal_delta_eband(); } @@ -923,6 +853,43 @@ void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish"); + // 6) calculate the local occupation number matrix and energy correction in + // DFT+U + if (PARAM.inp.dft_plus_u) + { + // only old DFT+U method should calculated energy correction in esolver, + // new DFT+U method will calculate energy in calculating Hamiltonian + if (PARAM.inp.dft_plus_u == 2) + { + if (GlobalC::dftu.omc != 2) + { + const std::vector>& tmp_dm + = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); + this->dftu_cal_occup_m(iter, tmp_dm); + } + GlobalC::dftu.cal_energy_correction(istep); + } + GlobalC::dftu.output(); + } + + // (7) for deepks, calculate delta_e +#ifdef __DEEPKS + if (PARAM.inp.deepks_scf) + { + const std::vector>& dm + = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); + + this->dpks_cal_e_delta_band(dm); + } +#endif + + // 8) for delta spin + if (PARAM.inp.sc_mag_switch) + { + SpinConstrain& sc = SpinConstrain::getScInstance(); + sc.cal_MW(iter, this->p_hamilt); + } + // call iter_finish() of ESolver_KS ESolver_KS::iter_finish(istep, iter); diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 5ebdd79548..b5f7c6093b 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -50,9 +50,7 @@ class ESolver_KS_LCAO : public ESolver_KS { virtual void iter_init(const int istep, const int iter) override; - virtual void hamilt2density(const int istep, - const int iter, - const double ethr) override; + virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; virtual void update_pot(const int istep, const int iter) override; diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 5c14fec3d0..0a7a88402b 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -118,10 +118,8 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell& this->pelec_td = dynamic_cast(this->pelec); } -void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, const double ethr) +void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int iter, const double ethr) { - pelec->charge->save_rho_before_sum_band(); - if (wf.init_wfc == "file") { if (istep >= 1) @@ -171,11 +169,23 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec_td, false); } } - // else - // { - // ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "HSolver has not been initialed!"); - // } + // symmetrize the charge density only for ground state + if (istep <= 1) + { + Symmetry_rho srho; + for (int is = 0; is < PARAM.inp.nspin; is++) + { + srho.begin(is, *(pelec->charge), pw_rho, GlobalC::ucell.symm); + } + } + + // (7) calculate delta energy + this->pelec->f_en.deband = this->pelec->cal_delta_eband(); +} + +void ESolver_KS_LCAO_TDDFT::iter_finish(const int istep, int& iter) +{ // print occupation of each band if (iter == 1 && istep <= 2) { @@ -201,32 +211,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons << std::endl; } - for (int ik = 0; ik < kv.get_nks(); ++ik) - { - this->pelec_td->print_band(ik, PARAM.inp.printe, iter); - } - - // using new charge density. - this->pelec->cal_energies(1); - - // symmetrize the charge density only for ground state - if (istep <= 1) - { - Symmetry_rho srho; - for (int is = 0; is < PARAM.inp.nspin; is++) - { - srho.begin(is, *(pelec->charge), pw_rho, GlobalC::ucell.symm); - } - } - - // (6) compute magnetization, only for spin==2 - GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, - this->pelec->charge->nxyz, - this->pelec->charge->rho, - pelec->nelec_spin.data()); - - // (7) calculate delta energy - this->pelec->f_en.deband = this->pelec->cal_delta_eband(); + ESolver_KS_LCAO, double>::iter_finish(istep, iter); } void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index 46a3ed084e..7589d77587 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -30,10 +30,12 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, doubl int td_htype = 1; protected: - virtual void hamilt2density(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; virtual void update_pot(const int istep, const int iter) override; + virtual void iter_finish(const int istep, int& iter) override; + virtual void after_scf(const int istep) override; void cal_edm_tddft(); diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp index 0d868dfc34..fc40a76eac 100644 --- a/source/module_esolver/esolver_ks_lcaopw.cpp +++ b/source/module_esolver/esolver_ks_lcaopw.cpp @@ -111,83 +111,55 @@ namespace ModuleESolver } template - void ESolver_KS_LIP::hamilt2density(const int istep, const int iter, const double ethr) + void ESolver_KS_LIP::hamilt2density_single(const int istep, const int iter, const double ethr) { - ModuleBase::TITLE("ESolver_KS_LIP", "hamilt2density"); - ModuleBase::timer::tick("ESolver_KS_LIP", "hamilt2density"); - + ModuleBase::TITLE("ESolver_KS_LIP", "hamilt2density_single"); + ModuleBase::timer::tick("ESolver_KS_LIP", "hamilt2density_single"); + + // reset energy + this->pelec->f_en.eband = 0.0; + this->pelec->f_en.demet = 0.0; + // choose if psi should be diag in subspace + // be careful that istep start from 0 and iter start from 1 + // if (iter == 1) + hsolver::DiagoIterAssist::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true; + hsolver::DiagoIterAssist::SCF_ITER = iter; + hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; + hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; + + // It is not a good choice to overload another solve function here, this will spoil the concept of + // multiple inheritance and polymorphism. But for now, we just do it in this way. + // In the future, there will be a series of class ESolver_KS_LCAO_PW, HSolver_LCAO_PW and so on. + std::weak_ptr> psig = this->p_wf_init->get_psig(); + + if (psig.expired()) { - // reset energy - this->pelec->f_en.eband = 0.0; - this->pelec->f_en.demet = 0.0; - // choose if psi should be diag in subspace - // be careful that istep start from 0 and iter start from 1 - // if (iter == 1) - hsolver::DiagoIterAssist::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true; - hsolver::DiagoIterAssist::SCF_ITER = iter; - hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; - - // It is not a good choice to overload another solve function here, this will spoil the concept of - // multiple inheritance and polymorphism. But for now, we just do it in this way. - // In the future, there will be a series of class ESolver_KS_LCAO_PW, HSolver_LCAO_PW and so on. - std::weak_ptr> psig = this->p_wf_init->get_psig(); - - if (psig.expired()) - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW::hamilt2density", "psig lifetime is expired"); - } + ModuleBase::WARNING_QUIT("ESolver_KS_PW::hamilt2density_single", "psig lifetime is expired"); + } - hsolver::HSolverLIP hsolver_lip_obj(this->pw_wfc); - hsolver_lip_obj.solve(this->p_hamilt, - this->kspw_psi[0], - this->pelec, - psig.lock().get()[0], - false); + hsolver::HSolverLIP hsolver_lip_obj(this->pw_wfc); + hsolver_lip_obj.solve(this->p_hamilt, this->kspw_psi[0], this->pelec, psig.lock().get()[0], false); - if (PARAM.inp.out_bandgap) - { - if (!PARAM.globalv.two_fermi) - { - this->pelec->cal_bandgap(); - } - else - { - this->pelec->cal_bandgap_updw(); - } - } - } - // add exx #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) { + if (GlobalC::exx_info.info_global.cal_exx) + { this->pelec->set_exx(this->exx_lip->get_exx_energy()); // Peize Lin add 2019-03-09 -} + } #endif - // calculate the delta_harris energy - // according to new charge density. - // mohan add 2009-01-23 - this->pelec->cal_energies(1); - Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::ucell.symm); } - // compute magnetization, only for LSDA(spin==2) - GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, - this->pelec->charge->nxyz, - this->pelec->charge->rho, - this->pelec->nelec_spin.data()); - // deband is calculated from "output" charge density calculated // in sum_band // need 'rho(out)' and 'vr (v_h(in) and v_xc(in))' this->pelec->f_en.deband = this->pelec->cal_delta_eband(); - ModuleBase::timer::tick("ESolver_KS_LIP", "hamilt2density"); + ModuleBase::timer::tick("ESolver_KS_LIP", "hamilt2density_single"); } template diff --git a/source/module_esolver/esolver_ks_lcaopw.h b/source/module_esolver/esolver_ks_lcaopw.h index c133791dcb..afb758848f 100644 --- a/source/module_esolver/esolver_ks_lcaopw.h +++ b/source/module_esolver/esolver_ks_lcaopw.h @@ -21,7 +21,7 @@ namespace ModuleESolver ~ESolver_KS_LIP(); /// All the other interfaces except this one are the same as ESolver_KS_PW. - virtual void hamilt2density(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; void before_all_runners(const Input_para& inp, UnitCell& cell) override; void after_all_runners()override; diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 91bc8c2618..51eb2c35a0 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -321,78 +321,49 @@ void ESolver_KS_PW::iter_init(const int istep, const int iter) // mohan move harris functional to here, 2012-06-05 // use 'rho(in)' and 'v_h and v_xc'(in) this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(); - - //(2) save change density as previous charge, - // prepared fox mixing. - if (GlobalV::MY_STOGROUP == 0) - { - this->pelec->charge->save_rho_before_sum_band(); - } } // Temporary, it should be replaced by hsolver later. template -void ESolver_KS_PW::hamilt2density(const int istep, const int iter, const double ethr) +void ESolver_KS_PW::hamilt2density_single(const int istep, const int iter, const double ethr) { - ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density"); - - { - // reset energy - this->pelec->f_en.eband = 0.0; - this->pelec->f_en.demet = 0.0; - // choose if psi should be diag in subspace - // be careful that istep start from 0 and iter start from 1 - // if (iter == 1) - hsolver::DiagoIterAssist::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true; - - hsolver::DiagoIterAssist::SCF_ITER = iter; - hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; - - hsolver::HSolverPW hsolver_pw_obj(this->pw_wfc, - &this->wf, - - PARAM.inp.calculation, - PARAM.inp.basis_type, - PARAM.inp.ks_solver, - PARAM.inp.use_paw, - PARAM.globalv.use_uspp, - PARAM.inp.nspin, - - hsolver::DiagoIterAssist::SCF_ITER, - hsolver::DiagoIterAssist::PW_DIAG_NMAX, - hsolver::DiagoIterAssist::PW_DIAG_THR, - - hsolver::DiagoIterAssist::need_subspace, - this->init_psi); - - hsolver_pw_obj.solve(this->p_hamilt, - this->kspw_psi[0], - this->pelec, - this->pelec->ekb.c, - GlobalV::RANK_IN_POOL, - GlobalV::NPROC_IN_POOL, - false); - - this->init_psi = true; - - if (PARAM.inp.out_bandgap) - { - if (!PARAM.globalv.two_fermi) - { - this->pelec->cal_bandgap(); - } - else - { - this->pelec->cal_bandgap_updw(); - } - } - } + ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density_single"); + + // reset energy + this->pelec->f_en.eband = 0.0; + this->pelec->f_en.demet = 0.0; + // choose if psi should be diag in subspace + // be careful that istep start from 0 and iter start from 1 + // if (iter == 1) + hsolver::DiagoIterAssist::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true; + + hsolver::DiagoIterAssist::SCF_ITER = iter; + hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; + hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; - // calculate the delta_harris energy - // according to new charge density. - // mohan add 2009-01-23 - this->pelec->cal_energies(1); + hsolver::HSolverPW hsolver_pw_obj(this->pw_wfc, + &this->wf, + PARAM.inp.calculation, + PARAM.inp.basis_type, + PARAM.inp.ks_solver, + PARAM.inp.use_paw, + PARAM.globalv.use_uspp, + PARAM.inp.nspin, + hsolver::DiagoIterAssist::SCF_ITER, + hsolver::DiagoIterAssist::PW_DIAG_NMAX, + hsolver::DiagoIterAssist::PW_DIAG_THR, + hsolver::DiagoIterAssist::need_subspace, + this->init_psi); + + hsolver_pw_obj.solve(this->p_hamilt, + this->kspw_psi[0], + this->pelec, + this->pelec->ekb.c, + GlobalV::RANK_IN_POOL, + GlobalV::NPROC_IN_POOL, + false); + + this->init_psi = true; Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) @@ -400,18 +371,12 @@ void ESolver_KS_PW::hamilt2density(const int istep, const int iter, c srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::ucell.symm); } - // compute magnetization, only for LSDA(spin==2) - GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, - this->pelec->charge->nxyz, - this->pelec->charge->rho, - this->pelec->nelec_spin.data()); - // deband is calculated from "output" charge density calculated // in sum_band // need 'rho(out)' and 'vr (v_h(in) and v_xc(in))' this->pelec->f_en.deband = this->pelec->cal_delta_eband(); - ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density"); + ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density_single"); } // Temporary, it should be rewritten with Hamilt class. diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index dc586a67a5..d27b1198c5 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -31,7 +31,7 @@ class ESolver_KS_PW : public ESolver_KS void cal_stress(ModuleBase::matrix& stress) override; - virtual void hamilt2density(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; virtual void hamilt2estates(const double ethr) override; diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index 023f800443..a0e93c325c 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -167,7 +167,7 @@ void ESolver_SDFT_PW::after_scf(const int istep) } template -void ESolver_SDFT_PW::hamilt2density(int istep, int iter, double ethr) +void ESolver_SDFT_PW::hamilt2density_single(int istep, int iter, double ethr) { // reset energy this->pelec->f_en.eband = 0.0; @@ -381,7 +381,7 @@ void ESolver_SDFT_PW::nscf() this->before_scf(istep); - this->hamilt2density(istep, iter, diag_thr); + this->hamilt2density_single(istep, iter, diag_thr); this->pelec->cal_energies(2); diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h index 63e613befc..9b00f33d78 100644 --- a/source/module_esolver/esolver_sdft_pw.h +++ b/source/module_esolver/esolver_sdft_pw.h @@ -35,7 +35,7 @@ class ESolver_SDFT_PW : public ESolver_KS_PW protected: virtual void before_scf(const int istep) override; - virtual void hamilt2density(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; virtual void nscf() override;