LCAO refactor step 4

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