diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 74a6e7bb17..3f76995e7a 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -253,7 +253,6 @@ OBJS_ESOLVER=esolver.o\ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\ esolver_ks_lcao_tddft.o\ dpks_cal_e_delta_band.o\ - dftu_cal_occup_m.o\ set_matrix_grid.o\ lcao_before_scf.o\ lcao_gets.o\ diff --git a/source/module_esolver/CMakeLists.txt b/source/module_esolver/CMakeLists.txt index 77795d82c0..17fb40401e 100644 --- a/source/module_esolver/CMakeLists.txt +++ b/source/module_esolver/CMakeLists.txt @@ -20,7 +20,6 @@ if(ENABLE_LCAO) esolver_ks_lcao_tddft.cpp dpks_cal_e_delta_band.cpp set_matrix_grid.cpp - dftu_cal_occup_m.cpp lcao_before_scf.cpp lcao_gets.cpp lcao_others.cpp diff --git a/source/module_esolver/dftu_cal_occup_m.cpp b/source/module_esolver/dftu_cal_occup_m.cpp deleted file mode 100644 index 5296305ec6..0000000000 --- a/source/module_esolver/dftu_cal_occup_m.cpp +++ /dev/null @@ -1,50 +0,0 @@ -#include "esolver_ks_lcao.h" -#include "module_hamilt_lcao/module_dftu/dftu.h" - -namespace ModuleESolver -{ - - using namespace std; - - //! dftu occupation matrix for gamma only using dm(double) - template <> - void ESolver_KS_LCAO::dftu_cal_occup_m( - const int& iter, - const vector>& dm)const - { - GlobalC::dftu.cal_occup_m_gamma( - iter, - dm, - this->p_chgmix->get_mixing_beta(), - this->p_hamilt); - } - - //! dftu occupation matrix for multiple k-points using dm(complex) - template <> - void ESolver_KS_LCAO, double>::dftu_cal_occup_m( - const int& iter, - const vector>>& dm)const - { - GlobalC::dftu.cal_occup_m_k( - iter, - dm, - this->kv, - this->p_chgmix->get_mixing_beta(), - this->p_hamilt); - } - - //! dftu occupation matrix - template <> - void ESolver_KS_LCAO, complex>::dftu_cal_occup_m( - const int& iter, - const vector>>& dm)const - { - GlobalC::dftu.cal_occup_m_k( - iter, - dm, - this->kv, - this->p_chgmix->get_mixing_beta(), - this->p_hamilt); - } - -} diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index cc0d56fc75..3d953fc4c5 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -22,19 +22,7 @@ namespace ModuleESolver { class ESolver_FP : public ESolver { - public: - - ModulePW::PW_Basis* pw_rho; - - /** - * @brief same as pw_rho for ncpp. Here 'd' stands for 'dense' - * dense grid for for uspp, used for ultrasoft augmented charge density. - * charge density and potential are defined on dense grids, - * but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi> - */ - ModulePW::PW_Basis* pw_rhod; - ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class - + public: //! Constructor ESolver_FP(); @@ -44,6 +32,10 @@ namespace ModuleESolver //! Initialize of the first-principels energy solver virtual void before_all_runners(const Input_para& inp, UnitCell& cell) override; + protected: + //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. + virtual void after_scf(const int istep); + virtual void init_after_vc(const Input_para& inp, UnitCell& cell); // liuyu add 2023-03-09 //! Electronic states @@ -58,9 +50,16 @@ namespace ModuleESolver //! K points in Brillouin zone K_Vectors kv; - protected: - //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. - virtual void after_scf(const int istep); + ModulePW::PW_Basis* pw_rho; + + /** + * @brief same as pw_rho for ncpp. Here 'd' stands for 'dense' + * dense grid for for uspp, used for ultrasoft augmented charge density. + * charge density and potential are defined on dense grids, + * but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi> + */ + ModulePW::PW_Basis* pw_rhod; + ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class //! Charge extrapolation Charge_Extra CE; diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 237cd7b36a..b088ebab02 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -5,6 +5,7 @@ #include "module_base/tool_title.h" #include "module_elecstate/module_dm/cal_dm_psi.h" #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" +#include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_io/berryphase.h" #include "module_io/cube_io.h" #include "module_io/dos_nao.h" @@ -879,7 +880,7 @@ void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) { const std::vector>& tmp_dm = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); - this->dftu_cal_occup_m(iter, tmp_dm); + ModuleDFTU::dftu_cal_occup_m(iter, tmp_dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt); } GlobalC::dftu.cal_energy_correction(istep); } diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 4acb416e63..7b2aa8cf71 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -97,10 +97,6 @@ class ESolver_KS_LCAO : public ESolver_KS { #endif private: - // tmp interfaces before sub-modules are refactored - void dftu_cal_occup_m(const int& iter, - const std::vector>& dm) const; - #ifdef __DEEPKS void dpks_cal_e_delta_band(const std::vector>& dm) const; diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index c77e857d76..5dac415a72 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -18,6 +18,15 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, doubl void before_all_runners(const Input_para& inp, UnitCell& cell) override; + protected: + 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; + //! wave functions of last time step psi::Psi>* psi_laststep = nullptr; @@ -31,16 +40,6 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, doubl elecstate::ElecStateLCAO_TDDFT* pelec_td = nullptr; int td_htype = 1; - - protected: - - 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; }; } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_ks_lcaopw.h b/source/module_esolver/esolver_ks_lcaopw.h index afb758848f..68fc0a905b 100644 --- a/source/module_esolver/esolver_ks_lcaopw.h +++ b/source/module_esolver/esolver_ks_lcaopw.h @@ -20,9 +20,6 @@ namespace ModuleESolver ~ESolver_KS_LIP(); - /// All the other interfaces except this one are the same as ESolver_KS_PW. - 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; @@ -30,6 +27,9 @@ namespace ModuleESolver virtual void iter_init(const int istep, const int iter) override; virtual void iter_finish(const int istep, int& iter) override; + /// All the other interfaces except this one are the same as ESolver_KS_PW. + virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; + virtual void allocate_hamilt() override; virtual void deallocate_hamilt() override; diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index f9476bafb8..160ee63c87 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -23,16 +23,12 @@ class ESolver_KS_PW : public ESolver_KS void before_all_runners(const Input_para& inp, UnitCell& cell) override; - void init_after_vc(const Input_para& inp, UnitCell& cell) override; - double cal_energy() override; void cal_force(ModuleBase::matrix& force) override; void cal_stress(ModuleBase::matrix& stress) override; - virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; - void after_all_runners() override; protected: @@ -48,6 +44,10 @@ class ESolver_KS_PW : public ESolver_KS virtual void others(const int istep) override; + virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; + + void init_after_vc(const Input_para& inp, UnitCell& cell) override; + // temporary, this will be removed in the future; // Init Global class void Init_GlobalC(const Input_para& inp, UnitCell& ucell, pseudopot_cell_vnl& ppcell); diff --git a/source/module_hamilt_lcao/module_dftu/dftu.cpp b/source/module_hamilt_lcao/module_dftu/dftu.cpp index 5bae255122..54c546c879 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu.cpp @@ -419,4 +419,25 @@ const hamilt::HContainer* DFTU::get_dmr(int ispin) const } } +//! dftu occupation matrix for gamma only using dm(double) +template <> +void dftu_cal_occup_m(const int iter, + const std::vector>& dm, + const K_Vectors& kv, + const double& mixing_beta, + hamilt::Hamilt* p_ham) +{ + GlobalC::dftu.cal_occup_m_gamma(iter, dm, mixing_beta, p_ham); +} + +//! dftu occupation matrix for multiple k-points using dm(complex) +template <> +void dftu_cal_occup_m(const int iter, + const std::vector>>& dm, + const K_Vectors& kv, + const double& mixing_beta, + hamilt::Hamilt>* p_ham) +{ + GlobalC::dftu.cal_occup_m_k(iter, dm, kv, mixing_beta, p_ham); +} } // namespace ModuleDFTU diff --git a/source/module_hamilt_lcao/module_dftu/dftu.h b/source/module_hamilt_lcao/module_dftu/dftu.h index c9b8eb34f9..fbd7e6dcfd 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.h +++ b/source/module_hamilt_lcao/module_dftu/dftu.h @@ -265,6 +265,13 @@ class DFTU const elecstate::DensityMatrix, double>* dm_in_dftu_cd = nullptr; }; +template +void dftu_cal_occup_m(const int iter, + const std::vector>& dm, + const K_Vectors& kv, + const double& mixing_beta, + hamilt::Hamilt* p_ham); + } // namespace ModuleDFTU namespace GlobalC