Refactor: move dftu_cal_occup_m to dftu

Author: YuLiu98

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/dftu_cal_occup_m.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"
@@ -873,7 +874,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.cpp

@@ -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.h

@@ -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.cpp

@@ -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