diff --git a/source/module_elecstate/elecstate_print.cpp b/source/module_elecstate/elecstate_print.cpp index 90a8598996..92c406a6b4 100644 --- a/source/module_elecstate/elecstate_print.cpp +++ b/source/module_elecstate/elecstate_print.cpp @@ -421,11 +421,6 @@ void ElecState::print_etot(const bool converged, {FmtTable::Align::LEFT, FmtTable::Align::CENTER}); table << titles << energies_Ry << energies_eV; GlobalV::ofs_running << table.str() << std::endl; - if (iter_in == 1) // pengfei Li added 2015-1-31 - { - this->f_en.etot_old = this->f_en.etot; - } - this->f_en.etot_delta = this->f_en.etot - this->f_en.etot_old; if (GlobalV::OUT_LEVEL == "ie" || GlobalV::OUT_LEVEL == "m") // xiaohui add 'm' option, 2015-09-16 { std::vector mag; @@ -462,7 +457,6 @@ void ElecState::print_etot(const bool converged, duration, 6); } - this->f_en.etot_old = this->f_en.etot; return; } diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 0fc82019e5..6525a4d32d 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -563,12 +563,6 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) } this->print_iter(iter, drho, dkin, duration, diag_ethr); - // add a energy threshold for SCF convergence - if (this->conv_elec == 0) // only check when density is not converged - { - this->conv_elec = ( iter != 1 && std::abs(this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < this->scf_ene_thr ); - } - // 12) Json, need to be moved to somewhere else #ifdef __RAPIDJSON // add Json of scf mag @@ -584,11 +578,7 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) if (this->conv_elec) { this->niter = iter; - bool stop = this->do_after_converge(iter); - if (stop) - { - break; - } + break; } // notice for restart @@ -615,6 +605,28 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) return; }; +template +void ESolver_KS::iter_finish(const int iter) +{ + // 1 means Harris-Foulkes functional + // 2 means Kohn-Sham functional + this->pelec->cal_energies(2); + + if (iter == 1) + { + this->pelec->f_en.etot_old = this->pelec->f_en.etot; + } + this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old; + this->pelec->f_en.etot_old = this->pelec->f_en.etot; + + // add a energy threshold for SCF convergence + if (this->conv_elec == 0) // only check when density is not converged + { + this->conv_elec + = (iter != 1 && std::abs(this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < this->scf_ene_thr); + } +} + //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. template void ESolver_KS::after_scf(const int istep) diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 2f42600144..558ab6c8e2 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -67,18 +67,15 @@ class ESolver_KS : public ESolver_FP virtual void iter_init(const int istep, const int iter) {}; //! Something to do after hamilt2density function in each iter loop. - virtual void iter_finish(const int iter) {}; + virtual void iter_finish(const int iter); - //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. + //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. virtual void after_scf(const int istep); //! It should be replaced by a function in Hamilt Class virtual void update_pot(const int istep, const int iter) {}; - //! choose strategy when charge density convergence achieved - virtual bool do_after_converge(int& iter){return true;} - - protected: + protected: // Print the headline on the screen: // ITER ETOT(eV) EDIFF(eV) DRHO TIME(s) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 7340a43e23..fdaed025b6 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -916,13 +916,15 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) //! 3) output exx matrix //! 4) output charge density and density matrix //! 5) cal_MW? (why put it here?) -//! 6) calculate the total energy? //------------------------------------------------------------------------------ template void ESolver_KS_LCAO::iter_finish(int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish"); + // call iter_finish() of ESolver_KS + ESolver_KS::iter_finish(iter); + // 1) mix density matrix if mixing_restart + mixing_dmr + not first // mixing_restart at every iter if (GlobalV::MIXING_RESTART > 0 && this->p_chgmix->mixing_restart_count > 0 && GlobalV::MIXING_DMR) @@ -980,16 +982,30 @@ void ESolver_KS_LCAO::iter_finish(int iter) GlobalC::restart.save_disk("Eexx", 0, 1, &this->pelec->f_en.exx); } } -#endif - // 4) output charge density and density matrix - bool print = false; - if (this->out_freq_elec && iter % this->out_freq_elec == 0) + if (GlobalC::exx_info.info_global.cal_exx && this->conv_elec) { - print = true; + if (GlobalC::exx_info.info_ri.real_number) + { + this->conv_elec = this->exd->exx_after_converge( + *this->p_hamilt, + *dynamic_cast*>(this->pelec)->get_DM(), + this->kv, + iter); + } + else + { + this->conv_elec = this->exc->exx_after_converge( + *this->p_hamilt, + *dynamic_cast*>(this->pelec)->get_DM(), + this->kv, + iter); + } } +#endif - if (print) + // 4) output charge density and density matrix + if (this->out_freq_elec && iter % this->out_freq_elec == 0) { for (int is = 0; is < GlobalV::NSPIN; is++) { @@ -1055,8 +1071,11 @@ void ESolver_KS_LCAO::iter_finish(int iter) sc.cal_MW(iter, this->p_hamilt); } - // 6) calculate the total energy. - this->pelec->cal_energies(2); + // 6) use the converged occupation matrix for next MD/Relax SCF calculation + if (GlobalV::dft_plus_u && this->conv_elec) + { + GlobalC::dftu.initialed_locale = true; + } } //------------------------------------------------------------------------------ @@ -1235,48 +1254,6 @@ void ESolver_KS_LCAO::after_scf(const int istep) } } -//------------------------------------------------------------------------------ -//! the 15th function of ESolver_KS_LCAO: do_after_converge -//! mohan add 2024-05-11 -//------------------------------------------------------------------------------ -template -bool ESolver_KS_LCAO::do_after_converge(int& iter) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "do_after_converge"); - - if (GlobalV::dft_plus_u) - { - // use the converged occupation matrix for next MD/Relax SCF calculation - GlobalC::dftu.initialed_locale = true; - } - // FIXME: for developer who want to test restarting DeePKS with same Descriptor/PDM in last MD step - // RUN: " GlobalC::ld.set_init_pdm(true); " can skip the calculation of PDM in the next iter_init - -#ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) - { - if (GlobalC::exx_info.info_ri.real_number) { - return this->exd->exx_after_converge( - *this->p_hamilt, - *dynamic_cast*>(this->pelec) - ->get_DM(), - this->kv, - iter); - } - else { - return this->exc->exx_after_converge( - *this->p_hamilt, - *dynamic_cast*>(this->pelec) - ->get_DM(), - this->kv, - iter); - } - } -#endif // __EXX - - return true; -} - //------------------------------------------------------------------------------ //! the 20th,21th,22th functions of ESolver_KS_LCAO //! mohan add 2024-05-11 diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 7e0791bb1e..c5c6d38047 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -60,8 +60,6 @@ class ESolver_KS_LCAO : public ESolver_KS { virtual void after_scf(const int istep) override; - virtual bool do_after_converge(int& iter) override; - virtual void others(const int istep) override; // we will get rid of this class soon, don't use it, mohan 2024-03-28 diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp index abb62cec87..b5542246b8 100644 --- a/source/module_esolver/esolver_ks_lcaopw.cpp +++ b/source/module_esolver/esolver_ks_lcaopw.cpp @@ -211,11 +211,13 @@ namespace ModuleESolver ModuleBase::timer::tick("ESolver_KS_LIP", "hamilt2density"); } -#ifdef __EXX template - bool ESolver_KS_LIP::do_after_converge(int& iter) + void ESolver_KS_LIP::iter_finish(const int iter) { - if (GlobalC::exx_info.info_global.cal_exx) + ESolver_KS_PW::iter_finish(iter); + +#ifdef __EXX + if (GlobalC::exx_info.info_global.cal_exx && this->conv_elec) { // no separate_loop case if (!GlobalC::exx_info.info_global.separate_loop) @@ -226,24 +228,24 @@ namespace ModuleESolver // in first scf loop, exx updated once in beginning, // in second scf loop, exx updated every iter - if (this->two_level_step) { - return true; - } else + if (!this->two_level_step) { // update exx and redo scf XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); iter = 0; std::cout << " Entering 2nd SCF, where EXX is updated" << std::endl; this->two_level_step++; - return false; + this->conv_elec = false; } } // has separate_loop case // exx converged or get max exx steps else if (this->two_level_step == GlobalC::exx_info.info_global.hybrid_step - || (iter == 1 && this->two_level_step != 0)) { - return true; - } else + || (iter == 1 && this->two_level_step != 0)) + { + this->conv_elec = true; + } + else { // update exx and redo scf if (this->two_level_step == 0) @@ -252,23 +254,24 @@ namespace ModuleESolver } std::cout << " Updating EXX " << std::flush; - timeval t_start; gettimeofday(&t_start, nullptr); + timeval t_start; + gettimeofday(&t_start, nullptr); this->exx_lip->cal_exx(); iter = 0; this->two_level_step++; - timeval t_end; gettimeofday(&t_end, nullptr); - std::cout << "and rerun SCF\t" - << std::setprecision(3) << std::setiosflags(std::ios::scientific) - << (double)(t_end.tv_sec - t_start.tv_sec) + (double)(t_end.tv_usec - t_start.tv_usec) / 1000000.0 - << std::defaultfloat << " (s)" << std::endl; - return false; + timeval t_end; + gettimeofday(&t_end, nullptr); + std::cout << "and rerun SCF\t" << std::setprecision(3) << std::setiosflags(std::ios::scientific) + << (double)(t_end.tv_sec - t_start.tv_sec) + + (double)(t_end.tv_usec - t_start.tv_usec) / 1000000.0 + << std::defaultfloat << " (s)" << std::endl; + this->conv_elec = false; } } - else { return true; } - } #endif + } template void ESolver_KS_LIP::after_all_runners() diff --git a/source/module_esolver/esolver_ks_lcaopw.h b/source/module_esolver/esolver_ks_lcaopw.h index 79bbd7cc23..dabe4822da 100644 --- a/source/module_esolver/esolver_ks_lcaopw.h +++ b/source/module_esolver/esolver_ks_lcaopw.h @@ -24,20 +24,19 @@ namespace ModuleESolver virtual void hamilt2density(const int istep, const int iter, const double ethr) override; void before_all_runners(const Input_para& inp, UnitCell& cell) override; - void iter_init(const int istep, const int iter) override; void after_all_runners()override; protected: - - virtual void allocate_hsolver() override; - virtual void deallocate_hsolver() override; - virtual void allocate_hamilt() override; - virtual void deallocate_hamilt() override; + virtual void iter_init(const int istep, const int iter) override; + virtual void iter_finish(const int iter) override; + virtual void allocate_hsolver() override; + virtual void deallocate_hsolver() override; + virtual void allocate_hamilt() override; + virtual void deallocate_hamilt() override; #ifdef __EXX std::unique_ptr> exx_lip; int two_level_step = 0; - virtual bool do_after_converge(int& iter) override; #endif }; diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index c9215eecb1..7717424ada 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -455,6 +455,9 @@ void ESolver_KS_PW::update_pot(const int istep, const int iter) template void ESolver_KS_PW::iter_finish(const int iter) { + // call iter_finish() of ESolver_KS + ESolver_KS::iter_finish(iter); + // liuyu 2023-10-24 // D in uspp need vloc, thus needs update when veff updated // calculate the effective coefficient matrix for non-local pseudopotential @@ -465,18 +468,7 @@ void ESolver_KS_PW::iter_finish(const int iter) GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, GlobalC::ucell); } - // 1 means Harris-Foulkes functional - // 2 means Kohn-Sham functional - const int energy_type = 2; - this->pelec->cal_energies(2); - - bool print = false; if (this->out_freq_elec && iter % this->out_freq_elec == 0) - { - print = true; - } - - if (print == true) { if (PARAM.inp.out_chg > 0) { diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index 42b4ba10e2..bbae24d243 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -139,7 +139,8 @@ void ESolver_SDFT_PW::before_scf(const int istep) void ESolver_SDFT_PW::iter_finish(int iter) { - this->pelec->cal_energies(2); + // call iter_finish() of ESolver_KS + ESolver_KS>::iter_finish(iter); } void ESolver_SDFT_PW::after_scf(const int istep)