diff --git a/source/Makefile.Objects b/source/Makefile.Objects index b05e69f8f9..e612261c6e 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -250,6 +250,7 @@ OBJS_ELECSTAT=elecstate.o\ read_pseudo.o\ cal_wfc.o\ setup_estate_pw.o\ + update_pot.o\ OBJS_ELECSTAT_LCAO=elecstate_lcao.o\ elecstate_lcao_cal_tau.o\ diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index a4e8126ad0..d9544209d3 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -23,6 +23,7 @@ #include "source_io/json_output/init_info.h" #include "source_io/json_output/output_info.h" +#include "source_estate/update_pot.h" // mohan add 20251016 namespace ModuleESolver { @@ -374,7 +375,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i if (this->scf_ene_thr > 0.0) { // calculate energy of output charge density - this->update_pot(ucell, istep, iter, conv_esolver); + elecstate::update_pot(ucell, this->pelec, this->chr, conv_esolver); 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; @@ -432,7 +433,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i #endif // 4) Update potentials (should be done every SF iter) - this->update_pot(ucell, istep, iter, conv_esolver); + elecstate::update_pot(ucell, this->pelec, this->chr, conv_esolver); // 5) calculate energies // 1 means Harris-Foulkes functional diff --git a/source/source_esolver/esolver_ks.h b/source/source_esolver/esolver_ks.h index d54794dd89..3f271dc4cc 100644 --- a/source/source_esolver/esolver_ks.h +++ b/source/source_esolver/esolver_ks.h @@ -46,9 +46,6 @@ class ESolver_KS : public ESolver_FP //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override; - //! It should be replaced by a function in Hamilt Class - virtual void update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver){}; - //! Hamiltonian hamilt::Hamilt* p_hamilt = nullptr; diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index b7601da811..c7171d0de7 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -446,22 +446,6 @@ void ESolver_KS_LCAO::hamilt2rho_single(UnitCell& ucell, int istep, int this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); } -template -void ESolver_KS_LCAO::update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) -{ - ModuleBase::TITLE("ESolver_KS_LCAO", "update_pot"); - - if (!conv_esolver) - { - elecstate::cal_ux(ucell); - this->pelec->pot->update_from_charge(&this->chr, &ucell); - this->pelec->f_en.descf = this->pelec->cal_delta_escf(); - } - else - { - this->pelec->cal_converged(); - } -} template void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) diff --git a/source/source_esolver/esolver_ks_lcao.h b/source/source_esolver/esolver_ks_lcao.h index 13359f7c8e..67ca7a0482 100644 --- a/source/source_esolver/esolver_ks_lcao.h +++ b/source/source_esolver/esolver_ks_lcao.h @@ -52,8 +52,6 @@ class ESolver_KS_LCAO : public ESolver_KS virtual void hamilt2rho_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override; - virtual void update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) override; - virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override; diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index c68ff0955b..3072579164 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -337,15 +337,19 @@ void ESolver_KS_LCAO_TDDFT::iter_finish( } ESolver_KS_LCAO, TR>::iter_finish(ucell, istep, iter, conv_esolver); + + this->save2(ucell, istep, iter, conv_esolver); + } template -void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, +void ESolver_KS_LCAO_TDDFT::save2(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) { // Calculate new potential according to new Charge Density +/* if (!conv_esolver) { elecstate::cal_ux(ucell); @@ -356,6 +360,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, { this->pelec->cal_converged(); } +*/ const int nloc = this->pv.nloc; const int ncol_nbands = this->pv.ncol_bands; diff --git a/source/source_esolver/esolver_ks_lcao_tddft.h b/source/source_esolver/esolver_ks_lcao_tddft.h index fc0a62e41d..75bd287591 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.h +++ b/source/source_esolver/esolver_ks_lcao_tddft.h @@ -64,7 +64,8 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, TR> virtual void hamilt2rho_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override; - virtual void update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) override; + // mohan change update_pot to save2, 2025-10-17 + void save2(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver); virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index d2e2612070..a8e5d553b7 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -28,6 +28,7 @@ #include "source_estate/setup_estate_pw.h" // mohan add 20251005 #include "source_io/ctrl_output_pw.h" // mohan add 20250927 #include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006 +#include "source_estate/update_pot.h" // mohan add 20251016 namespace ModuleESolver { @@ -268,24 +269,6 @@ void ESolver_KS_PW::hamilt2rho_single(UnitCell& ucell, const int iste ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2rho_single"); } -// Temporary, it should be rewritten with Hamilt class. -template -void ESolver_KS_PW::update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) -{ - if (!conv_esolver) - { - elecstate::cal_ux(ucell); - this->pelec->pot->update_from_charge(&this->chr, &ucell); - this->pelec->f_en.descf = this->pelec->cal_delta_escf(); -#ifdef __MPI - MPI_Bcast(&(this->pelec->f_en.descf), 1, MPI_DOUBLE, 0, BP_WORLD); -#endif - } - else - { - this->pelec->cal_converged(); - } -} template void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) @@ -343,8 +326,8 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int << std::endl; exx_helper.op_exx->first_iter = false; XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); - update_pot(ucell, istep, iter, conv_esolver); - exx_helper.iter_inc(); + elecstate::update_pot(ucell, this->pelec, this->chr, conv_esolver); + exx_helper.iter_inc(); } } } diff --git a/source/source_esolver/esolver_ks_pw.h b/source/source_esolver/esolver_ks_pw.h index acdd7083ee..09fbefc175 100644 --- a/source/source_esolver/esolver_ks_pw.h +++ b/source/source_esolver/esolver_ks_pw.h @@ -41,8 +41,6 @@ class ESolver_KS_PW : public ESolver_KS virtual void iter_init(UnitCell& ucell, const int istep, const int iter) override; - virtual void update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver) override; - virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override; diff --git a/source/source_esolver/rho_restart.cpp b/source/source_esolver/rho_restart.cpp index bd7497f34b..3a8c4b40b5 100644 --- a/source/source_esolver/rho_restart.cpp +++ b/source/source_esolver/rho_restart.cpp @@ -1,5 +1,6 @@ #include "rho_restart.h" +/* void ModuleESolver::rho_restart(const Input_para& inp, const UnitCell& ucell, const elecstate::ElecState& elec, @@ -145,3 +146,4 @@ void ModuleESolver::rho_restart(const Input_para& inp, #endif } +*/ diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index a8380ade4e..8a54110e55 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -39,6 +39,7 @@ list(APPEND objects read_pseudo.cpp cal_wfc.cpp setup_estate_pw.cpp + update_pot.cpp ) if(ENABLE_LCAO) diff --git a/source/source_estate/update_pot.cpp b/source/source_estate/update_pot.cpp new file mode 100644 index 0000000000..8a49f66e0c --- /dev/null +++ b/source/source_estate/update_pot.cpp @@ -0,0 +1,23 @@ +#include "source_estate/update_pot.h" +#include "source_estate/cal_ux.h" + +void elecstate::update_pot(UnitCell& ucell, // unitcell + elecstate::ElecState* &pelec, // pointer of electrons + const Charge &chr, + const bool conv_esolver + ) // charge density +{ + if (!conv_esolver) + { + elecstate::cal_ux(ucell); + pelec->pot->update_from_charge(&chr, &ucell); + pelec->f_en.descf = pelec->cal_delta_escf(); +#ifdef __MPI + MPI_Bcast(&(pelec->f_en.descf), 1, MPI_DOUBLE, 0, BP_WORLD); +#endif + } + else + { + pelec->cal_converged(); + } +} diff --git a/source/source_estate/update_pot.h b/source/source_estate/update_pot.h new file mode 100644 index 0000000000..cce4c42428 --- /dev/null +++ b/source/source_estate/update_pot.h @@ -0,0 +1,17 @@ +#ifndef UPDATE_POT_H +#define UPDATE_POT_H + +#include "source_cell/unitcell.h" +#include "source_estate/elecstate.h" + +namespace elecstate +{ + +void update_pot(UnitCell& ucell, // unitcell + elecstate::ElecState* &pelec, // pointer of electrons + const Charge &chr, + const bool conv_esolver); // charge density +} + + +#endif