LCAO refactor step 4 (#3921)

* update timer, add std

* update the format of esolver_fp.cpp

* add description in esolver_fp.h

* update a few esolver files

* update esolver description for LCAO

* update some formats of esolver codes

* update esolver_ks_pw.cpp formats

* update formats of esolver_ks_pw.cpp and esolver_ks_pw.h

* update esolver_ks_lcao_tddft.cpp formats

* update format of esolver_sdft_pw.cpp

* update the format of esolver_sdft_pw_tool.cpp

* keep formating esolver_ks_pw.cpp

* formating esolver_of_interface.cpp

* change GlobalC::ucell to ucell in esolver_ks.cpp

* remove some GlobalC::ucell in esolver_sdft_pw.cpp

* refactor the code before getting rid of RA in esolver_lcao

* refactor before getting rid of LOWF

* refactor before getting rid of LCAO_hamilt.h and LCAO_matrix.h

* refactor wavefunc_in_pw

* refactor density matrix

* refactor the format cal_dm_psi.cpp

* format forces.cpp

* refactor esolver_of_tool.cpp

* change member function beforescf in Esolver to before_scf

* change afterscf to after_scf

* change updatepot to update_pot

* change eachiterinit to iter_init, change eachiterfinish to iter_finish

* refactor esolvers, change member function names of most esolvers

* reformat esolver.h

* update tests for esolvers

* add TITLE in esolver_ks_lcao

* update esolver_ks_lcao

* update esolver_lcao

* update timer::tick in esolver_lcao

* try to delete LCAO_Matrix in LCAO_Hamilt, and try to delete Parallel_Orbitals in Force_k

* fix the compiling issue with LCAO_hamilt.hpp

* try to divide the FORCE_k.cpp into several small files

* divide FORCE_k into foverlap_k.cpp ftvnl_dphi_k.cpp fvl_dphi_k.cpp fvnl_dbeta_k.cpp four files

* get rid of UHM in FORCE_k.cpp

* cannot compile, but I have modified some files in order to get rid of Gint_k and Gint_Gamma in UHM

* keep updating, cannot run

* update write_Vxc, cannot run yet

* keep updating gint_gamma and gint_k

* update LCAO_matrix.cpp

* divide force files, and update Makefile.Objects

* update LCAO_hamilt.cpp

* delete genH pointer in UHM

* divide LCAO_hamilt.cpp into small codes, grid_init.cpp is the first one

* update grid_init in esolver_ks_lcao

* add a new namespace named sparse_format, most of the functions that originally belong to LCAO_hamilt should be moved to sparse_format

* cannot find the mismatch of DFTU

* fix the DFTU error

* update

* enable the test_memory function again by setting calcalculation parameter in INPUT file

* update memory record functions

* add sparse_format_u

* keep refactoring LCAO_hamilt

* rename sparse format files, which are originally defined in LCAO_hamilt.h

* add spar_u.h and spar_u.cpp

* update sparse matrix

* add spar_exx

* update LCAO_hamilt.cpp

* tear down LCAO_hamilt.h and .cpp, DONE.

* keep updating spar

* fix the bugs after deleting LCAO_hamilt

* continue

* fix some errors when compiling

* fix compiling errors in DOS

* fix errors without uhm

* remove directly use of LM in spar_hsr.cpp

* update spa_ files

* keep cleaning the mass left by UHM

* finnaly, all files can be compiled without uhm

* previous commit has problems, now it has been fixed.

* fix makefile bugs

* fix exx compiling errors

* fix exx problems

* update exx

* fix exx

* in some sense, the exx code is a disaster

* let's fix exx again

* fix undefination

* fix a bug in LCAO Refactor Step 4, I accidently forgot to add & for Grid_Driver& _grid in output_mat_sparse.h

---------

Co-authored-by: maki49 <1579492865@qq.com>

source/Makefile.Objects

@@ -501,9 +501,12 @@ OBJS_LCAO=DM_gamma.o\
       fvl_dphi_k.o\
       fvnl_dbeta_k.o\
       LCAO_gen_fixedH.o\
-      LCAO_hamilt.o\
       grid_init.o\
-      sparse_format.o\
+      spar_dh.o\
+      spar_exx.o\
+      spar_hsr.o\
+      spar_st.o\
+      spar_u.o\
       LCAO_matrix.o\
       LCAO_nnr.o\
       center2_orb-orb11.o\

source/module_elecstate/elecstate_lcao.h

@@ -21,7 +21,6 @@ class ElecStateLCAO : public ElecState
                   Local_Orbital_Charge* loc_in ,
                   Gint_Gamma* gint_gamma_in,  //mohan add 2024-04-01
                   Gint_k* gint_k_in, //mohan add 2024-04-01
-                  LCAO_Hamilt* uhm_in ,
                   Local_Orbital_wfc* lowf_in ,
                   ModulePW::PW_Basis* rhopw_in ,
                   ModulePW::PW_Basis_Big* bigpw_in )
@@ -77,7 +76,6 @@ class ElecStateLCAO : public ElecState
     Local_Orbital_Charge* loc = nullptr;
     Gint_Gamma* gint_gamma = nullptr; // mohan add 2024-04-01
     Gint_k* gint_k = nullptr; // mohan add 2024-04-01
-    //LCAO_Hamilt* uhm = nullptr; // mohan modify 2024-04-01
     Local_Orbital_wfc* lowf = nullptr;
     DensityMatrix<TK,double>* DM = nullptr;
 

source/module_esolver/esolver_ks_lcao.cpp

@@ -117,7 +117,6 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
             &(this->LOC),
             &(this->GG), // mohan add 2024-04-01
             &(this->GK), // mohan add 2024-04-01
-            &(this->uhm),
             &(this->LOWF),
             this->pw_rho,
             this->pw_big);
@@ -131,7 +130,7 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
     //------------------init Basis_lcao----------------------
 
     //! pass Hamilt-pointer to Operator
-    this->gen_h.LM = this->uhm.LM = &this->LM;
+    this->gen_h.LM = &this->LM;
 
     //! pass basis-pointer to EState and Psi
     this->LOC.ParaV = this->LOWF.ParaV = this->LM.ParaV = &(this->orb_con.ParaV);
@@ -259,7 +258,6 @@ void ESolver_KS_LCAO<TK, TR>::init_after_vc(Input& inp, UnitCell& ucell)
 				&(this->LOC),
                 &(this->GG), // mohan add 2024-04-01
                 &(this->GK), // mohan add 2024-04-01
-				&(this->uhm),
 				&(this->LOWF),
 				this->pw_rho,
 				this->pw_big);
@@ -426,7 +424,7 @@ void ESolver_KS_LCAO<TK, TR>::post_process(void)
     {
         ModuleIO::write_proj_band_lcao(
             this->psi,
-            this->uhm,
+            this->LM,
             this->pelec,
             this->kv,
             GlobalC::ucell,
@@ -437,7 +435,8 @@ void ESolver_KS_LCAO<TK, TR>::post_process(void)
     {
 		ModuleIO::out_dos_nao(
 				this->psi,
-				this->uhm,
+                this->LM,
+                this->orb_con.ParaV,
 				this->pelec->ekb,
 				this->pelec->wg,
 				INPUT.dos_edelta_ev,
@@ -1204,10 +1203,10 @@ ModuleIO::Output_Mat_Sparse<TK> ESolver_KS_LCAO<TK, TR>::create_Output_Mat_Spars
 			istep,
 			this->pelec->pot->get_effective_v(),
 			*this->LOWF.ParaV,
-			this->uhm,
-            this->gen_h, // mohan add 2024-04-02
+            this->gen_h, // mohan add 2024-04-06
             this->GK, // mohan add 2024-04-01
 			this->LM,
+            GlobalC::GridD, // mohan add 2024-04-06
 			this->kv,
 			this->p_hamilt);
 }

source/module_esolver/esolver_ks_lcao.h

@@ -5,7 +5,6 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h"
-#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h"
 // for grid integration
 #include "module_hamilt_lcao/module_gint/gint_gamma.h"
 #include "module_hamilt_lcao/module_gint/gint_k.h"
@@ -76,9 +75,6 @@ namespace ModuleESolver
         // we will get rid of this class soon, don't use it, mohan 2024-03-28
         Local_Orbital_Charge LOC;
 
-        // we will get rid of this class soon, don't use it, mohan 2024-03-28
-        LCAO_Hamilt uhm;
-
         LCAO_gen_fixedH gen_h; // mohan add 2024-04-02
 
         // used for k-dependent grid integration.
@@ -106,6 +102,7 @@ namespace ModuleESolver
         //--------------common for all calculation, not only scf-------------
         // set matrix and grid integral
         void set_matrix_grid(Record_adj& ra);
+
         void beforesolver(const int istep);
         //----------------------------------------------------------------------
 
@@ -127,11 +124,15 @@ namespace ModuleESolver
         std::shared_ptr<Exx_LRI<double>> exx_lri_double = nullptr;
         std::shared_ptr<Exx_LRI<std::complex<double>>> exx_lri_complex = nullptr;
 #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;
+
         void dpks_cal_projected_DM(const elecstate::DensityMatrix<TK, double>* dm) const;
 #endif
 

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -494,7 +494,9 @@ void ESolver_KS_LCAO<std::complex<double>, double>::get_S(void)
         dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, double>*>(this->p_hamilt->ops)->contributeHR();
     }
 
-    ModuleIO::output_S_R(this->uhm, this->p_hamilt, "SR.csr");
+    ModuleIO::output_SR(orb_con.ParaV, this->LM, GlobalC::GridD, this->p_hamilt, "SR.csr");
+
+    return;
 }
 
 
@@ -525,7 +527,10 @@ void ESolver_KS_LCAO<std::complex<double>, std::complex<double>>::get_S(void)
         dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>*>(this->p_hamilt->ops)
             ->contributeHR();
     }
-    ModuleIO::output_S_R(this->uhm, this->p_hamilt, "SR.csr");
+
+    ModuleIO::output_SR(orb_con.ParaV, this->LM, GlobalC::GridD, this->p_hamilt, "SR.csr");
+
+    return;
 }
 
 

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -100,7 +100,7 @@ void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell)
     //------------------init Hamilt_lcao----------------------
 
     // pass Hamilt-pointer to Operator
-    this->gen_h.LM = this->uhm.LM = &this->LM;
+    this->gen_h.LM = &this->LM;
     // pass basis-pointer to EState and Psi
     this->LOC.ParaV = this->LOWF.ParaV = this->LM.ParaV;
 

source/module_esolver/esolver_ks_lcao_tddft.h

@@ -4,7 +4,6 @@
 #include "esolver_ks_lcao.h"
 #include "module_basis/module_ao/ORB_control.h"
 #include "module_psi/psi.h"
-#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h"

source/module_esolver/esolver_ks_lcao_tmpfunc.cpp

@@ -4,67 +4,111 @@
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #endif
+
 namespace ModuleESolver
 {
+
+    using namespace std;
+
+    //! dftu occupation matrix for gamma only using dm(double)
     template <>
-    void ESolver_KS_LCAO<double, double>::dftu_cal_occup_m(const int& iter, const std::vector<std::vector<double>>& dm)const
+    void ESolver_KS_LCAO<double, double>::dftu_cal_occup_m(
+         const int& iter, 
+         const vector<vector<double>>& dm)const
     {
-        GlobalC::dftu.cal_occup_m_gamma(iter, dm, this->p_chgmix->get_mixing_beta());
+		GlobalC::dftu.cal_occup_m_gamma(
+				iter, 
+				dm, 
+				this->p_chgmix->get_mixing_beta());
     }
 
+    //! dftu occupation matrix for multiple k-points using dm(complex)
     template <>
-    void ESolver_KS_LCAO<std::complex<double>, double>::dftu_cal_occup_m(const int& iter, const std::vector<std::vector<std::complex<double>>>& dm)const
+    void ESolver_KS_LCAO<complex<double>, double>::dftu_cal_occup_m(
+       const int& iter, 
+       const vector<vector<complex<double>>>& dm)const
     {
-        GlobalC::dftu.cal_occup_m_k(iter, dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt);
+        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<std::complex<double>, std::complex<double>>::dftu_cal_occup_m(const int& iter, const std::vector<std::vector<std::complex<double>>>& dm)const
+    void ESolver_KS_LCAO<complex<double>, complex<double>>::dftu_cal_occup_m(
+       const int& iter, 
+       const vector<vector<complex<double>>>& dm)const
     {
-        GlobalC::dftu.cal_occup_m_k(iter, dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt);
+        GlobalC::dftu.cal_occup_m_k(
+                 iter, 
+                 dm, 
+                 this->kv, 
+                 this->p_chgmix->get_mixing_beta(), 
+                 this->p_hamilt);
     }
+
 #ifdef __DEEPKS
     template<>
-    void ESolver_KS_LCAO<double, double>::dpks_cal_e_delta_band(const std::vector<std::vector<double>>& dm)const
+    void ESolver_KS_LCAO<double, double>::dpks_cal_e_delta_band(
+       const vector<vector<double>>& dm)const
     {
         GlobalC::ld.cal_e_delta_band(dm);
     }
+
+
     template<>
-    void ESolver_KS_LCAO<std::complex<double>, double>::dpks_cal_e_delta_band(const std::vector<std::vector<std::complex<double>>>& dm)const
+    void ESolver_KS_LCAO<complex<double>, double>::dpks_cal_e_delta_band(
+        const vector<vector<complex<double>>>& dm)const
     {
         GlobalC::ld.cal_e_delta_band_k(dm, this->kv.nks);
     }
+
+
     template<>
-    void ESolver_KS_LCAO<std::complex<double>, std::complex<double>>::dpks_cal_e_delta_band(const std::vector<std::vector<std::complex<double>>>& dm)const
+    void ESolver_KS_LCAO<complex<double>, complex<double>>::dpks_cal_e_delta_band(
+        const vector<vector<complex<double>>>& dm)const
     {
         GlobalC::ld.cal_e_delta_band_k(dm, this->kv.nks);
     }
+
+
     template<>
-    void ESolver_KS_LCAO<double, double>::dpks_cal_projected_DM(const elecstate::DensityMatrix<double, double>* dm)const
+    void ESolver_KS_LCAO<double, double>::dpks_cal_projected_DM(
+        const elecstate::DensityMatrix<double, double>* dm)const
     {
-        GlobalC::ld.cal_projected_DM(dm, //this->LOC.dm_gamma,
+        GlobalC::ld.cal_projected_DM(dm,
             GlobalC::ucell,
             GlobalC::ORB,
             GlobalC::GridD);
     }
+
+
     template<>
-    void ESolver_KS_LCAO<std::complex<double>, double>::dpks_cal_projected_DM(const elecstate::DensityMatrix<std::complex<double>, double>* dm)const
+    void ESolver_KS_LCAO<complex<double>, double>::dpks_cal_projected_DM(
+          const elecstate::DensityMatrix<complex<double>, double>* dm)const
     {
-        GlobalC::ld.cal_projected_DM_k(dm, //this->LOC.dm_k,
+        GlobalC::ld.cal_projected_DM_k(dm,
             GlobalC::ucell,
             GlobalC::ORB,
             GlobalC::GridD,
             this->kv.nks,
             this->kv.kvec_d);
     }
+
+
     template<>
-    void ESolver_KS_LCAO<std::complex<double>, std::complex<double>>::dpks_cal_projected_DM(const elecstate::DensityMatrix<std::complex<double>, double>* dm)const
+    void ESolver_KS_LCAO<complex<double>, complex<double>>::dpks_cal_projected_DM(
+          const elecstate::DensityMatrix<complex<double>, double>* dm)const
     {
-        GlobalC::ld.cal_projected_DM_k(dm, //this->LOC.dm_k,
+        GlobalC::ld.cal_projected_DM_k(dm,
             GlobalC::ucell,
             GlobalC::ORB,
             GlobalC::GridD,
             this->kv.nks,
             this->kv.kvec_d);
     }
 #endif
-}
+}

source/module_hamilt_general/hamilt.h

@@ -21,10 +21,17 @@ class Hamilt
     virtual void updateHk(const int ik){return;}
 
     /// refresh status of Hamiltonian, for example, refresh H(R) and S(R) in LCAO case
-    virtual void refresh(){return;}
+    virtual void refresh(void){return;}
 
     /// core function: for solving eigenvalues of Hamiltonian with iterative method
-    virtual void hPsi(const T* psi_in, T* hpsi, const size_t size) const{return;}
+	virtual void hPsi(
+			const T* psi_in, 
+			T* hpsi, 
+			const size_t size) const
+	{
+		return;
+	}
+
     virtual void sPsi(const T* psi_in, // psi
                       T* spsi,         // spsi
                       const int nrow,  // dimension of spsi: nbands * nrow
@@ -35,21 +42,28 @@ class Hamilt
         syncmem_op()(this->ctx, this->ctx, spsi, psi_in, static_cast<size_t>(nbands * nrow));
     }
 
-    /// core function: return H(k) and S(k) matrixs for direct solving eigenvalues.
-    virtual void matrix(MatrixBlock<std::complex<double>> &hk_in, MatrixBlock<std::complex<double>> &sk_in){return;}
-    virtual void matrix(MatrixBlock<double> &hk_in, MatrixBlock<double> &sk_in){return;}
+	/// core function: return H(k) and S(k) matrixs for direct solving eigenvalues.
+	virtual void matrix(
+			MatrixBlock<std::complex<double>> &hk_in, 
+			MatrixBlock<std::complex<double>> &sk_in){return;}
+
+	virtual void matrix(
+			MatrixBlock<double> &hk_in, 
+			MatrixBlock<double> &sk_in){return;}
 
     std::string classname = "none";
 
     int non_first_scf=0;
 
     /// first node operator, add operations from each operators
     Operator<T, Device>* ops = nullptr;
+
 protected:
+
     Device* ctx = {};
     using syncmem_op = psi::memory::synchronize_memory_op<T, Device, Device>;
 };
 
 } // namespace hamilt
 
-#endif
+#endif

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -26,9 +26,12 @@ if(ENABLE_LCAO)
         fvl_dphi_k.cpp
         fvnl_dbeta_k.cpp
         LCAO_gen_fixedH.cpp
-        LCAO_hamilt.cpp
         grid_init.cpp
-        sparse_format.cpp
+        spar_dh.cpp
+        spar_exx.cpp
+        spar_hsr.cpp
+        spar_st.cpp
+        spar_u.cpp
         LCAO_matrix.cpp
         LCAO_nnr.cpp
         record_adj.cpp