Refactor: move dftu_cal_occup_m to dftu

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\

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

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;

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<TK, TR>::iter_finish(const int istep, int& iter)
             {
                 const std::vector<std::vector<TK>>& tmp_dm
                     = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(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);
         }

source/module_esolver/esolver_ks_lcao.h

@@ -97,10 +97,6 @@ class ESolver_KS_LCAO : public ESolver_KS<TK> {
 #endif
 
   private:
-    // tmp interfaces  before sub-modules are refactored
-    void dftu_cal_occup_m(const int& iter,
-                          const std::vector<std::vector<TK>>& dm) const;
-
 #ifdef __DEEPKS
     void dpks_cal_e_delta_band(const std::vector<std::vector<TK>>& dm) const;
 

source/module_esolver/esolver_ks_lcao_tddft.h

@@ -18,6 +18,15 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO<std::complex<double>, 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<std::complex<double>>* psi_laststep = nullptr;
 
@@ -31,16 +40,6 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO<std::complex<double>, 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

source/module_esolver/esolver_ks_lcaopw.h

@@ -20,16 +20,16 @@ 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;
 
     protected:
       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;
 

source/module_esolver/esolver_ks_pw.h

@@ -23,16 +23,12 @@ class ESolver_KS_PW : public ESolver_KS<T, Device>
 
     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<T, Device>
 
     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);

source/module_hamilt_lcao/module_dftu/dftu.cpp

@@ -419,4 +419,25 @@ const hamilt::HContainer<double>* 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<std::vector<double>>& dm,
+                      const K_Vectors& kv,
+                      const double& mixing_beta,
+                      hamilt::Hamilt<double>* 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<std::vector<std::complex<double>>>& dm,
+                      const K_Vectors& kv,
+                      const double& mixing_beta,
+                      hamilt::Hamilt<std::complex<double>>* p_ham)
+{
+    GlobalC::dftu.cal_occup_m_k(iter, dm, kv, mixing_beta, p_ham);
+}
 } // namespace ModuleDFTU

source/module_hamilt_lcao/module_dftu/dftu.h

@@ -265,6 +265,13 @@ class DFTU
     const elecstate::DensityMatrix<std::complex<double>, double>* dm_in_dftu_cd = nullptr;
 };
 
+template <typename T>
+void dftu_cal_occup_m(const int iter,
+                      const std::vector<std::vector<T>>& dm,
+                      const K_Vectors& kv,
+                      const double& mixing_beta,
+                      hamilt::Hamilt<T>* p_ham);
+
 } // namespace ModuleDFTU
 
 namespace GlobalC