diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index 82f74925303..aa1d627b825 100644
--- a/source/Makefile.Objects
+++ b/source/Makefile.Objects
@@ -502,6 +502,8 @@ OBJS_LCAO=DM_gamma.o\
       fvnl_dbeta_k.o\
       LCAO_gen_fixedH.o\
       LCAO_hamilt.o\
+      grid_init.o\
+      sparse_format.o\
       LCAO_matrix.o\
       LCAO_nnr.o\
       center2_orb-orb11.o\
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
index 01b98ce4923..c278f7e4c6a 100644
--- a/source/module_esolver/esolver_ks_lcao.cpp
+++ b/source/module_esolver/esolver_ks_lcao.cpp
@@ -117,7 +117,7 @@ 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->uhm),
             &(this->LOWF),
             this->pw_rho,
             this->pw_big);
@@ -131,7 +131,7 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
     //------------------init Basis_lcao----------------------
 
     //! pass Hamilt-pointer to Operator
-    this->UHM.genH.LM = this->UHM.LM = &this->LM;
+    this->gen_h.LM = this->uhm.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 +259,7 @@ 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->uhm),
 				&(this->LOWF),
 				this->pw_rho,
 				this->pw_big);
@@ -313,10 +313,10 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(ModuleBase::matrix& force)
 			GlobalV::TEST_STRESS,
 			this->LOC,
             this->orb_con.ParaV, 
-            this->LM,
 			this->pelec,
 			this->psi,
-			this->UHM,
+            this->LM,
+            this->gen_h, // mohan add 2024-04-02
             this->GG, // mohan add 2024-04-01
             this->GK, // mohan add 2024-04-01
 			force,
@@ -426,7 +426,7 @@ void ESolver_KS_LCAO<TK, TR>::post_process(void)
     {
         ModuleIO::write_proj_band_lcao(
             this->psi,
-            this->UHM,
+            this->uhm,
             this->pelec,
             this->kv,
             GlobalC::ucell,
@@ -437,7 +437,7 @@ void ESolver_KS_LCAO<TK, TR>::post_process(void)
     {
 		ModuleIO::out_dos_nao(
 				this->psi,
-				this->UHM,
+				this->uhm,
 				this->pelec->ekb,
 				this->pelec->wg,
 				INPUT.dos_edelta_ev,
@@ -1204,7 +1204,8 @@ 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->uhm,
+            this->gen_h, // mohan add 2024-04-02
             this->GK, // mohan add 2024-04-01
 			this->LM,
 			this->kv,
diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h
index c8c81ab203d..6e31bacd39c 100644
--- a/source/module_esolver/esolver_ks_lcao.h
+++ b/source/module_esolver/esolver_ks_lcao.h
@@ -77,7 +77,9 @@ namespace ModuleESolver
         Local_Orbital_Charge LOC;
 
         // we will get rid of this class soon, don't use it, mohan 2024-03-28
-        LCAO_Hamilt UHM;
+        LCAO_Hamilt uhm;
+
+        LCAO_gen_fixedH gen_h; // mohan add 2024-04-02
 
         // used for k-dependent grid integration.
         Gint_k GK;
diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp
index 09b6183ffb3..b01522210ee 100644
--- a/source/module_esolver/esolver_ks_lcao_elec.cpp
+++ b/source/module_esolver/esolver_ks_lcao_elec.cpp
@@ -23,6 +23,7 @@
 #include "module_io/dm_io.h"
 
 #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
 
 namespace ModuleESolver
 {
@@ -69,7 +70,7 @@ void ESolver_KS_LCAO<TK, TR>::set_matrix_grid(Record_adj& ra)
     // (2)For each atom, calculate the adjacent atoms in different cells
     // and allocate the space for H(R) and S(R).
     // If k point is used here, allocate HlocR after atom_arrange.
-    Parallel_Orbitals* pv = this->UHM.LM->ParaV;
+    Parallel_Orbitals* pv = this->LM.ParaV;
     ra.for_2d(*pv, GlobalV::GAMMA_ONLY_LOCAL);
     if (!GlobalV::GAMMA_ONLY_LOCAL)
     {
@@ -126,7 +127,7 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
     }
 
     // prepare grid in Gint
-	this->UHM.grid_prepare(
+	LCAO_domain::grid_prepare(
 			this->GridT, 
 			this->GG,
 			this->GK,
@@ -144,7 +145,7 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
         elecstate::DensityMatrix<TK, double>* DM = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
         this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr,
             GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK),
-            &(this->UHM.genH),
+            &(this->gen_h),
             &(this->LM),
             &(this->LOC),
             this->pelec->pot,
@@ -234,7 +235,7 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
     // since it depends on ionic positions
     if (GlobalV::deepks_setorb)
     {
-        const Parallel_Orbitals* pv = this->UHM.LM->ParaV;
+        const Parallel_Orbitals* pv = this->LM.ParaV;
         // build and save <psi(0)|alpha(R)> at beginning
         GlobalC::ld.build_psialpha(GlobalV::CAL_FORCE, GlobalC::ucell, GlobalC::ORB, GlobalC::GridD, GlobalC::UOT);
 
@@ -476,14 +477,15 @@ void ESolver_KS_LCAO<std::complex<double>, double>::get_S()
                          GlobalV::test_atom_input);
 
     this->RA.for_2d(this->orb_con.ParaV, GlobalV::GAMMA_ONLY_LOCAL);
-    this->UHM.genH.LM->ParaV = &this->orb_con.ParaV;
+    this->LM.ParaV = &this->orb_con.ParaV;
     if (this->p_hamilt == nullptr)
     {
-        this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(this->UHM.genH.LM, this->kv);
+        this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(&this->LM, this->kv);
         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_S_R(this->uhm, this->p_hamilt, "SR.csr");
 }
+
 template <>
 void ESolver_KS_LCAO<std::complex<double>, std::complex<double>>::get_S()
 {
@@ -503,15 +505,15 @@ void ESolver_KS_LCAO<std::complex<double>, std::complex<double>>::get_S()
                          GlobalV::test_atom_input);
 
     this->RA.for_2d(this->orb_con.ParaV, GlobalV::GAMMA_ONLY_LOCAL);
-    this->UHM.genH.LM->ParaV = &this->orb_con.ParaV;
+    this->LM.ParaV = &this->orb_con.ParaV;
     if (this->p_hamilt == nullptr)
     {
         this->p_hamilt
-            = new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(this->UHM.genH.LM, this->kv);
+            = new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(&this->LM, this->kv);
         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_S_R(this->uhm, this->p_hamilt, "SR.csr");
 }
 
 template <typename TK, typename TR>
diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp
index 182242eb6a3..f104a2454cc 100644
--- a/source/module_esolver/esolver_ks_lcao_tddft.cpp
+++ b/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->UHM.genH.LM = this->UHM.LM = &this->LM;
+    this->gen_h.LM = this->uhm.LM = &this->LM;
     // pass basis-pointer to EState and Psi
     this->LOC.ParaV = this->LOWF.ParaV = this->LM.ParaV;
 
@@ -455,8 +455,9 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep)
                         kv,
                         tmp_DM->get_paraV_pointer(),
                         this->RA,
-                        this->UHM);
-    }
+						this->LM, // mohan add 2024-04-02
+						this->gen_h); // mohan add 2024-02
+	}
     ESolver_KS_LCAO<std::complex<double>, double>::after_scf(istep);
 }
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt
index 61a21928f1f..2de30cca691 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt
@@ -27,6 +27,8 @@ if(ENABLE_LCAO)
         fvnl_dbeta_k.cpp
         LCAO_gen_fixedH.cpp
         LCAO_hamilt.cpp
+        grid_init.cpp
+        sparse_format.cpp
         LCAO_matrix.cpp
         LCAO_nnr.cpp
         record_adj.cpp
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
index 91f0e2665a1..4e98f79458b 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
@@ -31,10 +31,10 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                                           const bool istests,
 										  Local_Orbital_Charge &loc,
 										  Parallel_Orbitals &pv,
-										  LCAO_Matrix &lm,
                                           const elecstate::ElecState* pelec,
                                           const psi::Psi<T>* psi,
-										  LCAO_Hamilt& uhm,
+										  LCAO_Matrix& lm,
+                                          LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 										  Gint_Gamma &gint_gamma, // mohan add 2024-04-01
 										  Gint_k &gint_k, // mohan add 2024-04-01
 										  ModuleBase::matrix& fcs,
@@ -159,7 +159,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
 #else
                                      svl_dphi,
 #endif
-									 uhm,
+                                     gen_h, // mohan add 2024-04-02
 									 gint_gamma,
 									 gint_k,
 									 pv,
@@ -236,18 +236,19 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
         }
         if(GlobalV::dft_plus_u == 2)
         {
-            GlobalC::dftu.force_stress(pelec, *uhm.LM, force_dftu, stress_dftu, kv);
+            GlobalC::dftu.force_stress(pelec, lm, force_dftu, stress_dftu, kv);
         }
         else
         {
-            hamilt::DFTU<hamilt::OperatorLCAO<T, double>> tmp_dftu(uhm.LM,
-                                                                 kv.kvec_d,
-                                                                 nullptr,
-                                                                 nullptr,
-                                                                 GlobalC::ucell,
-                                                                 &GlobalC::GridD,
-                                                                 &GlobalC::dftu,
-                                                                 *(uhm.LM->ParaV));
+			hamilt::DFTU<hamilt::OperatorLCAO<T, double>> tmp_dftu(
+					&lm,
+					kv.kvec_d,
+					nullptr,
+					nullptr,
+					GlobalC::ucell,
+					&GlobalC::GridD,
+					&GlobalC::dftu,
+					*(lm.ParaV));
             tmp_dftu.cal_force_stress(isforce, isstress, force_dftu, stress_dftu);
         }
     }
@@ -743,7 +744,7 @@ void Force_Stress_LCAO<double>::integral_part(
 #else
     ModuleBase::matrix& svl_dphi,
 #endif
-	LCAO_Hamilt &uhm,
+    LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 	Gint_Gamma &gint_gamma, // mohan add 2024-04-01
 	Gint_k &gint_k, // mohan add 2024-04-01
 	Parallel_Orbitals &pv,
@@ -769,7 +770,7 @@ void Force_Stress_LCAO<double>::integral_part(
 #else
         svl_dphi,
 #endif
-        uhm,
+        gen_h,
         gint_gamma,
         lm);
     return;
@@ -797,7 +798,7 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(
 #else
     ModuleBase::matrix& svl_dphi,
 #endif
-	LCAO_Hamilt &uhm,
+    LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 	Gint_Gamma &gint_gamma,
 	Gint_k &gint_k,
 	Parallel_Orbitals &pv,
@@ -823,7 +824,7 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(
 #else
                      svl_dphi,
 #endif
-					 uhm,
+					 gen_h,
                      gint_k,
 					 pv,
 					 lm,
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
index 1604fd76cc8..f59bb821115 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
@@ -34,10 +34,10 @@ class Force_Stress_LCAO
         const bool istests,
 		Local_Orbital_Charge& loc,
 		Parallel_Orbitals &pv,
-		LCAO_Matrix &lm,
 		const elecstate::ElecState* pelec,
         const psi::Psi<T>* psi,
-		LCAO_Hamilt& uhm,
+		LCAO_Matrix &lm,
+        LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 		Gint_Gamma &gint_gamma, // mohan add 2024-04-01
 		Gint_k &gint_k, // mohan add 2024-04-01
         ModuleBase::matrix& fcs,
@@ -92,7 +92,7 @@ class Force_Stress_LCAO
 #else
         ModuleBase::matrix& svl_dphi,
 #endif
-		LCAO_Hamilt &uhm,
+        LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 		Gint_Gamma &gint_gamma,
 		Gint_k &gint_k,
 	    Parallel_Orbitals &pv,
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
index 0b5835446b1..f54a0ca85c9 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
@@ -38,7 +38,7 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce,
 #else
                                     ModuleBase::matrix& svl_dphi,
 #endif
-									LCAO_Hamilt &uhm,
+		                            LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
                                     Gint_Gamma &gint_gamma,
 									LCAO_Matrix &lm)
 {
@@ -54,7 +54,7 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce,
 
     // allocate DSloc_x, DSloc_y, DSloc_z
     // allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z
-    this->allocate_gamma(*this->ParaV, uhm.genH, lm);
+    this->allocate_gamma(*this->ParaV, gen_h, lm);
 
     // calculate the 'energy density matrix' here.
     this->cal_foverlap(isforce, isstress, psid, pelec, lm, foverlap, soverlap);
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h
index c8fe00c040c..d1c99b4dca2 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h
@@ -45,7 +45,7 @@ class Force_LCAO_gamma
 #else
                       ModuleBase::matrix& svl_dphi,
 #endif
-					  LCAO_Hamilt& uhm,
+                      LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
 					  Gint_Gamma &gint_gamma,
 					  LCAO_Matrix& lm);
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
index 98fd9379793..18a8dd4ad9b 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
@@ -50,7 +50,7 @@ void Force_LCAO_k::ftable_k(const bool isforce,
 #else
                             ModuleBase::matrix& svl_dphi,
 #endif
-                            LCAO_Hamilt &uhm,
+                            LCAO_gen_fixedH &gen_h,
                             Gint_k &gint_k,
                             Parallel_Orbitals &pv,
                             LCAO_Matrix &lm,
@@ -65,7 +65,7 @@ void Force_LCAO_k::ftable_k(const bool isforce,
 	this->allocate_k(
 			pv, 
 			lm, 
-            uhm.genH,
+            gen_h,
 			kv.nks, 
 			kv.kvec_d);
 
@@ -148,7 +148,7 @@ void Force_LCAO_k::ftable_k(const bool isforce,
             GlobalC::ld.check_v_delta_k(pv->nnr);
             for (int ik = 0; ik < kv.nks; ik++)
             {
-                uhm.LM->folding_fixedH(ik, kv.kvec_d);
+                LM->folding_fixedH(ik, kv.kvec_d);
             }
             GlobalC::ld.cal_e_delta_band_k(dm_k, kv.nks);
             std::ofstream ofs("E_delta_bands.dat");
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h
index 9fa60e2ad5e..34f09eadaf0 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h
@@ -42,7 +42,7 @@ class Force_LCAO_k : public Force_LCAO_gamma
 #else
                   ModuleBase::matrix& svl_dphi,
 #endif
-				  LCAO_Hamilt &uhm,
+                  LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
                   Gint_k &gint_k,
 				  Parallel_Orbitals &pv,
 				  LCAO_Matrix &lm,
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h
new file mode 100644
index 00000000000..25133c8b52b
--- /dev/null
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h
@@ -0,0 +1,24 @@
+#ifndef LCAO_DOMAIN_H
+#define LCAO_DOMAIN_H
+
+#include "module_cell/module_neighbor/sltk_atom_arrange.h"
+#include "module_cell/module_neighbor/sltk_grid_driver.h"
+#include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
+
+
+namespace LCAO_domain
+{
+
+//! prepare grid integration 
+void grid_prepare(
+		const Grid_Technique& gt, 
+        Gint_Gamma &gint_gamma,
+        Gint_k &gint_k,
+		const ModulePW::PW_Basis& rhopw, 
+		const ModulePW::PW_Basis_Big& bigpw);
+
+
+}
+
+#endif
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.h
index 24357e7a1df..913d514d93b 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.h
@@ -13,9 +13,6 @@
 
 class LCAO_gen_fixedH
 {
-    friend class Force_LCAO_gamma;
-    friend class Force_LCAO_k;
-    friend class LCAO_Hamilt;
 
   public:
     LCAO_Matrix* LM;
@@ -37,7 +34,6 @@ class LCAO_gen_fixedH
                       double dmax = 0.0);
 	// cal_syns : calculate asynchronous overlap matrix for Hefei-NAMD
 
-  private:
     // can used in gamma algorithm.
     void build_Nonlocal_beta_new(double* Hloc);
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp
index b1881849f48..998311df6c7 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp
@@ -16,6 +16,10 @@
 #include "LCAO_hamilt.hpp"
 #endif
 
+#include "sparse_format.h"
+
+using namespace sparse_format;
+
 LCAO_Hamilt::LCAO_Hamilt()
 {
 }
@@ -28,187 +32,6 @@ LCAO_Hamilt::~LCAO_Hamilt()
     }
 }
 
-//--------------------------------------------
-// prepare grid network for Gint(grid integral)
-//--------------------------------------------
-void LCAO_Hamilt::grid_prepare(
-		const Grid_Technique& gt, 
-        Gint_Gamma &gint_gamma,
-        Gint_k &gint_k,
-		const ModulePW::PW_Basis& rhopw, 
-		const ModulePW::PW_Basis_Big& bigpw)
-{
-    ModuleBase::TITLE("LCAO_Hamilt","grid_prepare");
-    ModuleBase::timer::tick("LCAO_Hamilt","grid_prepare");
-
-    if(GlobalV::GAMMA_ONLY_LOCAL)
-    {
-        gint_gamma.prep_grid(gt, bigpw.nbx, bigpw.nby, bigpw.nbzp, bigpw.nbzp_start,
-            rhopw.nxyz, bigpw.bx, bigpw.by, bigpw.bz, bigpw.bxyz, bigpw.nbxx,
-            rhopw.ny, rhopw.nplane, rhopw.startz_current);
-
-    }
-    else // multiple k-points
-    {
-        // cal the grid integration of 'Vl' matrix for l-points algorithms.
-        gint_k.prep_grid(gt, bigpw.nbx, bigpw.nby, bigpw.nbzp, bigpw.nbzp_start,
-            rhopw.nxyz, bigpw.bx, bigpw.by, bigpw.bz, bigpw.bxyz, bigpw.nbxx,
-            rhopw.ny, rhopw.nplane, rhopw.startz_current);
-    }
-
-    ModuleBase::timer::tick("LCAO_Hamilt","grid_prepare");
-    return;
-}
-
-void LCAO_Hamilt::set_R_range_sparse(LCAO_Matrix &lm)
-{
-    int R_minX = int(GlobalC::GridD.getD_minX());
-    int R_minY = int(GlobalC::GridD.getD_minY());
-    int R_minZ = int(GlobalC::GridD.getD_minZ());
-
-    int R_x = GlobalC::GridD.getCellX();
-    int R_y = GlobalC::GridD.getCellY();
-    int R_z = GlobalC::GridD.getCellZ();
-
-    for(int ix = 0; ix < R_x; ix++)
-    {
-        for(int iy = 0; iy < R_y; iy++)
-        {
-            for(int iz = 0; iz < R_z; iz++)
-            {
-                Abfs::Vector3_Order<int> temp_R(ix+R_minX, iy+R_minY, iz+R_minZ);
-                lm.all_R_coor.insert(temp_R);
-            }
-        }
-    }
-
-    return;
-}
-
-void LCAO_Hamilt::cal_STN_R_sparse_for_S(
-		const Parallel_Orbitals pv,
-        std::vector<double> &slocR,
-        std::vector<std::complex<double>> &slocR_soc,
-        std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>> &SR_sparse,
-        std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> &SR_soc_sparse,
-		const double &sparse_threshold)
-{
-    ModuleBase::TITLE("LCAO_Hamilt","cal_STN_R_sparse_for_S");
-
-    int index = 0;
-    ModuleBase::Vector3<double> dtau, tau1, tau2;
-    ModuleBase::Vector3<double> dtau1, dtau2, tau0;
-
-    double tmp=0.0;
-    std::complex<double> tmpc=complex<double>(0.0,0.0); 
-
-    for(int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1)
-    {
-        Atom* atom1 = &GlobalC::ucell.atoms[T1];
-        for(int I1 = 0; I1 < atom1->na; ++I1)
-        {
-            tau1 = atom1->tau[I1];
-            GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1);
-            Atom* atom1 = &GlobalC::ucell.atoms[T1];
-            const int start = GlobalC::ucell.itiaiw2iwt(T1,I1,0);
-
-            for(int ad = 0; ad < GlobalC::GridD.getAdjacentNum()+1; ++ad)
-            {
-                const int T2 = GlobalC::GridD.getType(ad);
-                const int I2 = GlobalC::GridD.getNatom(ad);
-                Atom* atom2 = &GlobalC::ucell.atoms[T2];
-
-                tau2 = GlobalC::GridD.getAdjacentTau(ad);
-                dtau = tau2 - tau1;
-                double distance = dtau.norm() * GlobalC::ucell.lat0;
-                double rcut = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ORB.Phi[T2].getRcut();
-
-                bool adj = false;
-
-				if(distance < rcut) 
-				{
-					adj = true;
-				}
-                else if(distance >= rcut)
-                {
-                    for(int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0)
-                    {
-                        const int T0 = GlobalC::GridD.getType(ad0);
-
-                        tau0 = GlobalC::GridD.getAdjacentTau(ad0);
-                        dtau1 = tau0 - tau1;
-                        dtau2 = tau0 - tau2;
-
-                        double distance1 = dtau1.norm() * GlobalC::ucell.lat0;
-                        double distance2 = dtau2.norm() * GlobalC::ucell.lat0;
-
-                        double rcut1 = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
-                        double rcut2 = GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
-
-                        if( distance1 < rcut1 && distance2 < rcut2 )
-                        {
-                            adj = true;
-                            break;
-                        }
-                    }
-                }
-
-                if(adj)
-                {
-                    const int start2 = GlobalC::ucell.itiaiw2iwt(T2,I2,0);
-
-					Abfs::Vector3_Order<int> dR(
-							GlobalC::GridD.getBox(ad).x, 
-							GlobalC::GridD.getBox(ad).y, 
-							GlobalC::GridD.getBox(ad).z);
-
-                    for(int ii=0; ii<atom1->nw*GlobalV::NPOL; ii++)
-                    {
-                        const int iw1_all = start + ii;
-                        const int mu = pv.global2local_row(iw1_all);
-
-						if(mu<0)
-						{
-							continue;
-						}
-
-                        for(int jj=0; jj<atom2->nw*GlobalV::NPOL; jj++)
-                        {
-                            int iw2_all = start2 + jj;
-                            const int nu = pv.global2local_col(iw2_all);
-
-							if(nu<0)
-							{
-								continue;
-							}
-
-                            if(GlobalV::NSPIN!=4)
-                            {
-                                tmp = slocR[index];
-                                if (std::abs(tmp) > sparse_threshold)
-                                {
-                                    SR_sparse[dR][iw1_all][iw2_all] = tmp;
-                                }
-                            }
-                            else
-                            {
-                                tmpc = slocR_soc[index];
-                                if(std::abs(tmpc) > sparse_threshold)
-                                {
-                                    SR_soc_sparse[dR][iw1_all][iw2_all] = tmpc;
-                                }
-                            }
-
-                            ++index;
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    return;
-}
 
 #include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
 void LCAO_Hamilt::cal_HContainer_sparse_d(
@@ -306,7 +129,7 @@ void LCAO_Hamilt::cal_HSR_sparse(
 {
     ModuleBase::TITLE("LCAO_Hamilt","cal_HSR_sparse");
 
-    set_R_range_sparse(*this->LM);
+    sparse_format::set_R_range(*this->LM);
 
     //cal_STN_R_sparse(current_spin, sparse_threshold);
     if(GlobalV::NSPIN!=4)
@@ -372,47 +195,6 @@ void LCAO_Hamilt::cal_HSR_sparse(
     clear_zero_elements(current_spin, sparse_threshold);
 }
 
-void LCAO_Hamilt::cal_dH_sparse(
- const int &current_spin, 
- const double &sparse_threshold,
- Gint_k &gint_k)
-{
-    ModuleBase::TITLE("LCAO_Hamilt","cal_dH_sparse");
-
-    set_R_range_sparse(*this->LM);
-
-    const int nnr = this->LM->ParaV->nnr;
-    this->LM->DHloc_fixedR_x = new double[nnr];
-    this->LM->DHloc_fixedR_y = new double[nnr];
-    this->LM->DHloc_fixedR_z = new double[nnr];
-
-    ModuleBase::GlobalFunc::ZEROS(this->LM->DHloc_fixedR_x, this->LM->ParaV->nloc);
-    ModuleBase::GlobalFunc::ZEROS(this->LM->DHloc_fixedR_y, this->LM->ParaV->nloc);
-    ModuleBase::GlobalFunc::ZEROS(this->LM->DHloc_fixedR_z, this->LM->ParaV->nloc);
-    // cal dT=<phi|kin|dphi> in LCAO
-    // cal T + VNL(P1) in LCAO basis
-    if(GlobalV::CAL_STRESS)
-	{
-        GlobalV::CAL_STRESS = false;
-        this->genH.build_ST_new('T', true, GlobalC::ucell, this->LM->Hloc_fixedR.data());
-        GlobalV::CAL_STRESS = true;
-    }
-    else
-    {
-        this->genH.build_ST_new('T', true, GlobalC::ucell, this->LM->Hloc_fixedR.data());
-    }
-    this->genH.build_Nonlocal_mu_new (this->LM->Hloc_fixed.data(), true);
-    
-    cal_dSTN_R_sparse(current_spin, sparse_threshold);
-
-    delete[] this->LM->DHloc_fixedR_x;
-    delete[] this->LM->DHloc_fixedR_y;
-    delete[] this->LM->DHloc_fixedR_z;
-
-    gint_k.cal_dvlocal_R_sparseMatrix(current_spin, sparse_threshold, this->LM);
-}
-
-
 void LCAO_Hamilt::cal_STN_R_sparse_for_T(const double &sparse_threshold)
 {
     ModuleBase::TITLE("LCAO_Hamilt","cal_STN_R_sparse_for_T");
@@ -526,7 +308,7 @@ void LCAO_Hamilt::cal_STN_R_sparse_for_T(const double &sparse_threshold)
 void LCAO_Hamilt::cal_SR_sparse(const double &sparse_threshold, hamilt::Hamilt<std::complex<double>>* p_ham)
 {
     ModuleBase::TITLE("LCAO_Hamilt","cal_SR_sparse");
-    set_R_range_sparse(*this->LM);
+    sparse_format::set_R_range(*this->LM);
     //cal_STN_R_sparse(current_spin, sparse_threshold);
     if(GlobalV::NSPIN!=4)
     {
@@ -542,15 +324,19 @@ void LCAO_Hamilt::cal_SR_sparse(const double &sparse_threshold, hamilt::Hamilt<s
     }
 }
 
-void LCAO_Hamilt::cal_TR_sparse(const double &sparse_threshold)
+void LCAO_Hamilt::cal_TR_sparse(
+		LCAO_gen_fixedH &gen_h,
+		const double &sparse_threshold)
 {
     ModuleBase::TITLE("LCAO_Hamilt","cal_TR_sparse");
     
     //need to rebuild T(R)
     this->LM->Hloc_fixedR.resize(this->LM->ParaV->nnr);
     this->LM->zeros_HSR('T');
-    this->genH.build_ST_new('T', 0, GlobalC::ucell, this->LM->Hloc_fixedR.data());
-    this->set_R_range_sparse(*this->LM);
+
+    gen_h.build_ST_new('T', 0, GlobalC::ucell, this->LM->Hloc_fixedR.data());
+
+    sparse_format::set_R_range(*this->LM);
     this->cal_STN_R_sparse_for_T(sparse_threshold);
 
     return;
@@ -872,117 +658,6 @@ void LCAO_Hamilt::clear_zero_elements(const int &current_spin, const double &spa
     }
 }
 
-void LCAO_Hamilt::cal_dSTN_R_sparse(const int &current_spin, const double &sparse_threshold)
-{
-    ModuleBase::TITLE("LCAO_Hamilt","cal_dSTN_R_sparse");
-
-    int index = 0;
-    ModuleBase::Vector3<double> dtau, tau1, tau2;
-    ModuleBase::Vector3<double> dtau1, dtau2, tau0;
-
-    double temp_value_double;
-    std::complex<double> temp_value_complex;
-
-    for(int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1)
-    {
-        Atom* atom1 = &GlobalC::ucell.atoms[T1];
-        for(int I1 = 0; I1 < atom1->na; ++I1)
-        {
-            tau1 = atom1->tau[I1];
-            GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1);
-            Atom* atom1 = &GlobalC::ucell.atoms[T1];
-            const int start = GlobalC::ucell.itiaiw2iwt(T1,I1,0);
-
-            for(int ad = 0; ad < GlobalC::GridD.getAdjacentNum()+1; ++ad)
-            {
-                const int T2 = GlobalC::GridD.getType(ad);
-                const int I2 = GlobalC::GridD.getNatom(ad);
-                Atom* atom2 = &GlobalC::ucell.atoms[T2];
-
-                tau2 = GlobalC::GridD.getAdjacentTau(ad);
-                dtau = tau2 - tau1;
-                double distance = dtau.norm() * GlobalC::ucell.lat0;
-                double rcut = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ORB.Phi[T2].getRcut();
-
-                bool adj = false;
-
-                if(distance < rcut) adj = true;
-                else if(distance >= rcut)
-                {
-                    for(int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0)
-                    {
-                        const int T0 = GlobalC::GridD.getType(ad0);
-
-                        tau0 = GlobalC::GridD.getAdjacentTau(ad0);
-                        dtau1 = tau0 - tau1;
-                        dtau2 = tau0 - tau2;
-
-                        double distance1 = dtau1.norm() * GlobalC::ucell.lat0;
-                        double distance2 = dtau2.norm() * GlobalC::ucell.lat0;
-
-                        double rcut1 = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
-                        double rcut2 = GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
-
-                        if( distance1 < rcut1 && distance2 < rcut2 )
-                        {
-                            adj = true;
-                            break;
-                        }
-                    }
-                }
-
-                if(adj)
-                {
-                    const int start2 = GlobalC::ucell.itiaiw2iwt(T2,I2,0);
-
-                    Abfs::Vector3_Order<int> dR(GlobalC::GridD.getBox(ad).x, GlobalC::GridD.getBox(ad).y, GlobalC::GridD.getBox(ad).z);
-
-                    for(int ii=0; ii<atom1->nw*GlobalV::NPOL; ii++)
-                    {
-                        const int iw1_all = start + ii;
-                        const int mu = this->LM->ParaV->global2local_row(iw1_all);
-
-                        if(mu<0)continue;
-
-                        for(int jj=0; jj<atom2->nw*GlobalV::NPOL; jj++)
-                        {
-                            int iw2_all = start2 + jj;
-                            const int nu = this->LM->ParaV->global2local_col(iw2_all);
-
-                            if(nu<0)continue;
-
-                            if(GlobalV::NSPIN!=4)
-                            {
-                                temp_value_double = this->LM->DHloc_fixedR_x[index];
-                                if (std::abs(temp_value_double) > sparse_threshold)
-                                {
-                                    this->LM->dHRx_sparse[current_spin][dR][iw1_all][iw2_all] = temp_value_double;
-                                }
-                                temp_value_double = this->LM->DHloc_fixedR_y[index];
-                                if (std::abs(temp_value_double) > sparse_threshold)
-                                {
-                                    this->LM->dHRy_sparse[current_spin][dR][iw1_all][iw2_all] = temp_value_double;
-                                }
-                                temp_value_double = this->LM->DHloc_fixedR_z[index];
-                                if (std::abs(temp_value_double) > sparse_threshold)
-                                {
-                                    this->LM->dHRz_sparse[current_spin][dR][iw1_all][iw2_all] = temp_value_double;
-                                }
-                            }
-                            else
-                            {
-                                ModuleBase::WARNING_QUIT("cal_dSTN_R_sparse","soc not supported!");
-                            }
-                            ++index;
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    return;
-}
 
 void LCAO_Hamilt::destroy_all_HSR_sparse(void)
 {
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h
index b57ebf43e7b..df6f5dfc165 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h
@@ -24,16 +24,6 @@ class LCAO_Hamilt
 
     ~LCAO_Hamilt();
 
-	void grid_prepare(
-			const Grid_Technique& gt, 
-			Gint_Gamma &gint_gamma,
-			Gint_k &gint_k,
-			const ModulePW::PW_Basis& rhopw, 
-			const ModulePW::PW_Basis_Big& bigpw);
-
-    // jingan add 2021-6-4
-    void set_R_range_sparse(LCAO_Matrix &lm);
-
     void cal_HContainer_sparse_d(const int &current_spin, 
         const double &sparse_threshold, 
         const hamilt::HContainer<double>& hR, 
@@ -44,16 +34,6 @@ class LCAO_Hamilt
         const hamilt::HContainer<std::complex<double>>& hR, 
         std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>>& target);
 
-    void cal_dSTN_R_sparse(const int &current_spin, const double &sparse_threshold);
-
-    // mohan update 2024-04-01
-	void cal_STN_R_sparse_for_S(
-			const Parallel_Orbitals pv,
-			std::vector<double> &slocR,
-			std::vector<std::complex<double>> &slocR_soc,
-			std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>> &SR_sparse,
-			std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> &SR_soc_sparse,
-			const double &sparse_threshold);
 
     void cal_STN_R_sparse_for_T(const double &sparse_threshold);
 
@@ -81,20 +61,14 @@ class LCAO_Hamilt
 
     void destroy_all_HSR_sparse(void);
 
-    void cal_TR_sparse(const double &sparse_threshold);
+    void cal_TR_sparse(
+			LCAO_gen_fixedH &gen_h,
+			const double &sparse_threshold);
 
     void destroy_TR_sparse(void);
 
-	void cal_dH_sparse(
-			const int &current_spin, 
-			const double &sparse_threshold,
-			Gint_k &gint_k); // mohan add 2024-04-01
-
     void destroy_dH_R_sparse(void);
 
-    // use overlap matrix to generate fixed Hamiltonian
-    LCAO_gen_fixedH genH;
-
     LCAO_Matrix* LM;
 };
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/grid_init.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/grid_init.cpp
new file mode 100644
index 00000000000..a385c14d035
--- /dev/null
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/grid_init.cpp
@@ -0,0 +1,63 @@
+#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
+#include "module_base/parallel_reduce.h"
+#include "module_base/timer.h"
+
+namespace LCAO_domain
+{
+
+//--------------------------------------------
+// prepare grid network for Gint(grid integral)
+//--------------------------------------------
+void grid_prepare(
+		const Grid_Technique& gt, 
+        Gint_Gamma &gint_gamma,
+        Gint_k &gint_k,
+		const ModulePW::PW_Basis& rhopw, 
+		const ModulePW::PW_Basis_Big& bigpw)
+{
+    ModuleBase::TITLE("LCAO_domain","grid_prepare");
+    ModuleBase::timer::tick("LCAO_domain","grid_prepare");
+
+    if(GlobalV::GAMMA_ONLY_LOCAL)
+    {
+		gint_gamma.prep_grid(
+				gt, 
+				bigpw.nbx, 
+				bigpw.nby, 
+				bigpw.nbzp, 
+				bigpw.nbzp_start,
+				rhopw.nxyz, 
+				bigpw.bx, 
+				bigpw.by, 
+				bigpw.bz, 
+				bigpw.bxyz, 
+				bigpw.nbxx,
+				rhopw.ny, 
+				rhopw.nplane, 
+				rhopw.startz_current);
+	}
+    else // multiple k-points
+    {
+        // cal the grid integration of 'Vl' matrix for l-points algorithms.
+		gint_k.prep_grid(
+				gt, 
+				bigpw.nbx, 
+				bigpw.nby, 
+				bigpw.nbzp, 
+				bigpw.nbzp_start,
+				rhopw.nxyz, 
+				bigpw.bx, 
+				bigpw.by, 
+				bigpw.bz, 
+				bigpw.bxyz, 
+				bigpw.nbxx,
+				rhopw.ny, 
+				rhopw.nplane, 
+				rhopw.startz_current);
+	}
+
+    ModuleBase::timer::tick("LCAO_domain","grid_prepare");
+    return;
+}
+
+}
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp
index b2490f18ab0..e6bfbc4c783 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp
@@ -121,7 +121,7 @@ void Veff<OperatorLCAO<TK, TR>>::contributeHR()
 
 // special case of gamma-only
 template<>
-void Veff<OperatorLCAO<double, double>>::contributeHR()
+void Veff<OperatorLCAO<double, double>>::contributeHR(void)
 {
     ModuleBase::TITLE("Veff", "contributeHR");
     ModuleBase::timer::tick("Veff", "contributeHR");
@@ -160,5 +160,4 @@ template class Veff<OperatorLCAO<double, double>>;
 template class Veff<OperatorLCAO<std::complex<double>, double>>;
 
 template class Veff<OperatorLCAO<std::complex<double>, std::complex<double>>>;
-
 }
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/sparse_format.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/sparse_format.cpp
new file mode 100644
index 00000000000..640c57f055e
--- /dev/null
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/sparse_format.cpp
@@ -0,0 +1,196 @@
+#include "sparse_format.h"
+
+void sparse_format::cal_dH(
+		LCAO_Matrix &lm,
+		LCAO_gen_fixedH &gen_h, 
+		const int &current_spin, 
+		const double &sparse_threshold,
+		Gint_k &gint_k)
+{
+    ModuleBase::TITLE("sparse_format","cal_dH");
+
+    sparse_format::set_R_range(lm);
+
+    const int nnr = lm.ParaV->nnr;
+    lm.DHloc_fixedR_x = new double[nnr];
+    lm.DHloc_fixedR_y = new double[nnr];
+    lm.DHloc_fixedR_z = new double[nnr];
+
+    ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixedR_x, lm.ParaV->nloc);
+    ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixedR_y, lm.ParaV->nloc);
+    ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixedR_z, lm.ParaV->nloc);
+    // cal dT=<phi|kin|dphi> in LCAO
+    // cal T + VNL(P1) in LCAO basis
+    if(GlobalV::CAL_STRESS)
+	{
+        GlobalV::CAL_STRESS = false;
+
+        gen_h.build_ST_new('T', true, GlobalC::ucell, lm.Hloc_fixedR.data());
+
+        GlobalV::CAL_STRESS = true;
+    }
+    else
+    {
+        gen_h.build_ST_new('T', true, GlobalC::ucell, lm.Hloc_fixedR.data());
+    }
+    gen_h.build_Nonlocal_mu_new (lm.Hloc_fixed.data(), true);
+    
+    sparse_format::cal_dSTN_R(lm, current_spin, sparse_threshold);
+
+    delete[] lm.DHloc_fixedR_x;
+    delete[] lm.DHloc_fixedR_y;
+    delete[] lm.DHloc_fixedR_z;
+
+    gint_k.cal_dvlocal_R_sparseMatrix(current_spin, sparse_threshold, &lm);
+}
+
+
+void sparse_format::set_R_range(LCAO_Matrix &lm)
+{
+    int R_minX = int(GlobalC::GridD.getD_minX());
+    int R_minY = int(GlobalC::GridD.getD_minY());
+    int R_minZ = int(GlobalC::GridD.getD_minZ());
+
+    int R_x = GlobalC::GridD.getCellX();
+    int R_y = GlobalC::GridD.getCellY();
+    int R_z = GlobalC::GridD.getCellZ();
+
+    for(int ix = 0; ix < R_x; ix++)
+    {
+        for(int iy = 0; iy < R_y; iy++)
+        {
+            for(int iz = 0; iz < R_z; iz++)
+            {
+                Abfs::Vector3_Order<int> temp_R(ix+R_minX, iy+R_minY, iz+R_minZ);
+                lm.all_R_coor.insert(temp_R);
+            }
+        }
+    }
+
+    return;
+}
+
+
+void sparse_format::cal_dSTN_R(
+		LCAO_Matrix &lm,
+		const int &current_spin, 
+		const double &sparse_threshold)
+{
+    ModuleBase::TITLE("sparse_format","cal_dSTN_R");
+
+    int index = 0;
+    ModuleBase::Vector3<double> dtau, tau1, tau2;
+    ModuleBase::Vector3<double> dtau1, dtau2, tau0;
+
+    double temp_value_double;
+    std::complex<double> temp_value_complex;
+
+    for(int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1)
+    {
+        Atom* atom1 = &GlobalC::ucell.atoms[T1];
+        for(int I1 = 0; I1 < atom1->na; ++I1)
+        {
+            tau1 = atom1->tau[I1];
+            GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1);
+            Atom* atom1 = &GlobalC::ucell.atoms[T1];
+            const int start = GlobalC::ucell.itiaiw2iwt(T1,I1,0);
+
+            for(int ad = 0; ad < GlobalC::GridD.getAdjacentNum()+1; ++ad)
+            {
+                const int T2 = GlobalC::GridD.getType(ad);
+                const int I2 = GlobalC::GridD.getNatom(ad);
+                Atom* atom2 = &GlobalC::ucell.atoms[T2];
+
+                tau2 = GlobalC::GridD.getAdjacentTau(ad);
+                dtau = tau2 - tau1;
+                double distance = dtau.norm() * GlobalC::ucell.lat0;
+                double rcut = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ORB.Phi[T2].getRcut();
+
+                bool adj = false;
+
+                if(distance < rcut) adj = true;
+                else if(distance >= rcut)
+                {
+                    for(int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0)
+                    {
+                        const int T0 = GlobalC::GridD.getType(ad0);
+
+                        tau0 = GlobalC::GridD.getAdjacentTau(ad0);
+                        dtau1 = tau0 - tau1;
+                        dtau2 = tau0 - tau2;
+
+                        double distance1 = dtau1.norm() * GlobalC::ucell.lat0;
+                        double distance2 = dtau2.norm() * GlobalC::ucell.lat0;
+
+                        double rcut1 = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
+                        double rcut2 = GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
+
+                        if( distance1 < rcut1 && distance2 < rcut2 )
+                        {
+                            adj = true;
+                            break;
+                        }
+                    }
+                }
+
+                if(adj)
+                {
+                    const int start2 = GlobalC::ucell.itiaiw2iwt(T2,I2,0);
+
+					Abfs::Vector3_Order<int> dR(
+							GlobalC::GridD.getBox(ad).x, 
+							GlobalC::GridD.getBox(ad).y, 
+							GlobalC::GridD.getBox(ad).z);
+
+                    for(int ii=0; ii<atom1->nw*GlobalV::NPOL; ii++)
+                    {
+                        const int iw1_all = start + ii;
+                        const int mu = lm.ParaV->global2local_row(iw1_all);
+
+						if(mu<0)
+						{
+							continue;
+						}
+
+                        for(int jj=0; jj<atom2->nw*GlobalV::NPOL; jj++)
+                        {
+                            int iw2_all = start2 + jj;
+                            const int nu = lm.ParaV->global2local_col(iw2_all);
+
+							if(nu<0)
+							{
+								continue;
+							}
+
+                            if(GlobalV::NSPIN!=4)
+                            {
+                                temp_value_double = lm.DHloc_fixedR_x[index];
+                                if (std::abs(temp_value_double) > sparse_threshold)
+                                {
+                                    lm.dHRx_sparse[current_spin][dR][iw1_all][iw2_all] = temp_value_double;
+                                }
+                                temp_value_double = lm.DHloc_fixedR_y[index];
+                                if (std::abs(temp_value_double) > sparse_threshold)
+                                {
+                                    lm.dHRy_sparse[current_spin][dR][iw1_all][iw2_all] = temp_value_double;
+                                }
+                                temp_value_double = lm.DHloc_fixedR_z[index];
+                                if (std::abs(temp_value_double) > sparse_threshold)
+                                {
+                                    lm.dHRz_sparse[current_spin][dR][iw1_all][iw2_all] = temp_value_double;
+                                }
+                            }
+                            else
+                            {
+                                ModuleBase::WARNING_QUIT("cal_dSTN_R","nspin=4 with SOC is not supported yet.");
+                            }
+                            ++index;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    return;
+}
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/sparse_format.h b/source/module_hamilt_lcao/hamilt_lcaodft/sparse_format.h
new file mode 100644
index 00000000000..1d7a09af780
--- /dev/null
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/sparse_format.h
@@ -0,0 +1,32 @@
+#ifndef SPARSE_FORMAT_H 
+#define SPARSE_FORMAT_H
+
+#include "module_cell/module_neighbor/sltk_atom_arrange.h"
+#include "module_cell/module_neighbor/sltk_grid_driver.h"
+#include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
+
+
+namespace sparse_format
+{
+
+void cal_dH(
+		LCAO_Matrix &lm,
+		LCAO_gen_fixedH &gen_h, 
+		const int &current_spin, 
+		const double &sparse_threshold,
+		Gint_k &gint_k);
+
+// be called by 'cal_dH_sparse'
+void set_R_range(LCAO_Matrix &lm);
+
+// be called by 'cal_dH_sparse' 
+void cal_dSTN_R(
+		LCAO_Matrix &lm,
+		const int &current_spin, 
+		const double &sparse_threshold);
+
+
+}
+
+#endif
diff --git a/source/module_io/cal_r_overlap_R.cpp b/source/module_io/cal_r_overlap_R.cpp
index 2743c99b6c6..9def5c977c7 100644
--- a/source/module_io/cal_r_overlap_R.cpp
+++ b/source/module_io/cal_r_overlap_R.cpp
@@ -319,13 +319,28 @@ void cal_r_overlap_R::out_rR(const int &istep)
                             int iL2 = iw2iL[orb_index_col];  
                             int im2 = iw2im[orb_index_col];
 
-                            ModuleBase::Vector3<double> r_distance = ( GlobalC::ucell.atoms[it2].tau[ia2] - GlobalC::ucell.atoms[it1].tau[ia1] + R_car ) * GlobalC::ucell.lat0;
+                            ModuleBase::Vector3<double> r_distance = 
+                            ( GlobalC::ucell.atoms[it2].tau[ia2] - 
+                            GlobalC::ucell.atoms[it1].tau[ia1] + R_car ) * GlobalC::ucell.lat0;
 
-                            double overlap_o = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( origin_point, r_distance, im1, im2 );
-                            double overlap_x = -1 * factor * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( origin_point, r_distance, im1, 1, im2 ); // m =  1
-                            double overlap_y = -1 * factor * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( origin_point, r_distance, im1, 2, im2 ); // m = -1
-                            double overlap_z =      factor * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( origin_point, r_distance, im1, 0, im2 ); // m =  0
-                            ModuleBase::Vector3<double> temp_prp = ModuleBase::Vector3<double>(overlap_x, overlap_y, overlap_z) + GlobalC::ucell.atoms[it1].tau[ia1] * GlobalC::ucell.lat0 * overlap_o;
+                            double overlap_o = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( 
+                            origin_point, r_distance, im1, im2 );
+
+                            double overlap_x = -1 * factor * 
+                            center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( 
+                            origin_point, r_distance, im1, 1, im2 ); // m =  1
+
+                            double overlap_y = -1 * factor * 
+                            center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( 
+                            origin_point, r_distance, im1, 2, im2 ); // m = -1
+
+                            double overlap_z = factor * 
+                            center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( 
+                            origin_point, r_distance, im1, 0, im2 ); // m =  0
+
+                            ModuleBase::Vector3<double> temp_prp = 
+                            ModuleBase::Vector3<double>(overlap_x, overlap_y, overlap_z) 
+                            + GlobalC::ucell.atoms[it1].tau[ia1] * GlobalC::ucell.lat0 * overlap_o;
 
                             if (std::abs(temp_prp.x) > sparse_threshold)
                             {
@@ -409,9 +424,13 @@ void cal_r_overlap_R::out_rR(const int &istep)
         std::ofstream out_r;
         std::stringstream ssr;
         if(GlobalV::CALCULATION == "md" && !GlobalV::out_app_flag)
-            ssr << GlobalV::global_matrix_dir << step << "_" << "data-rR-sparse.csr";
-        else
-            ssr << GlobalV::global_out_dir << "data-rR-sparse.csr";
+		{
+			ssr << GlobalV::global_matrix_dir << step << "_" << "data-rR-sparse.csr";
+		}
+		else
+		{
+			ssr << GlobalV::global_out_dir << "data-rR-sparse.csr";
+		}
 
         if (binary)
         {
@@ -478,10 +497,14 @@ void cal_r_overlap_R::out_rR_other(const int &istep, const std::set<Abfs::Vector
 
     std::ofstream out_r;
     std::stringstream ssr;
-    if(GlobalV::CALCULATION == "md" && !GlobalV::out_app_flag)
-        ssr << GlobalV::global_matrix_dir << step << "_" << "data-rR-sparse.csr";
-    else
-        ssr << GlobalV::global_out_dir << "data-rR-sparse.csr";
+	if(GlobalV::CALCULATION == "md" && !GlobalV::out_app_flag)
+	{
+		ssr << GlobalV::global_matrix_dir << step << "_" << "data-rR-sparse.csr";
+	}
+	else
+	{
+		ssr << GlobalV::global_out_dir << "data-rR-sparse.csr";
+	}
 
     if (GlobalV::DRANK == 0)
     {
@@ -525,7 +548,8 @@ void cal_r_overlap_R::out_rR_other(const int &istep, const std::set<Abfs::Vector
 
         ModuleBase::Vector3<double> R_car = ModuleBase::Vector3<double>(dRx, dRy, dRz) * GlobalC::ucell.latvec;
 
-        int ir, ic;
+        int ir=0; 
+        int ic=0;
         for(int iw1 = 0; iw1 < GlobalV::NLOCAL; iw1++)
         {
             ir = this->ParaV->global2local_row(iw1);
@@ -540,7 +564,8 @@ void cal_r_overlap_R::out_rR_other(const int &istep, const std::set<Abfs::Vector
                         int orb_index_col = iw2 / GlobalV::NPOL;
                         
                         // The off-diagonal term in SOC calculaiton is zero, and the two diagonal terms are the same
-                        int new_index = iw1 - GlobalV::NPOL*orb_index_row + (iw2 - GlobalV::NPOL*orb_index_col)*GlobalV::NPOL;
+                        int new_index = iw1 - GlobalV::NPOL*orb_index_row 
+                         + (iw2 - GlobalV::NPOL*orb_index_col)*GlobalV::NPOL;
                         
                         if(new_index == 0 || new_index == 3)
                         {
@@ -556,13 +581,29 @@ void cal_r_overlap_R::out_rR_other(const int &istep, const std::set<Abfs::Vector
                             int iL2 = iw2iL[orb_index_col];  
                             int im2 = iw2im[orb_index_col];
 
-                            ModuleBase::Vector3<double> r_distance = ( GlobalC::ucell.atoms[it2].tau[ia2] - GlobalC::ucell.atoms[it1].tau[ia1] + R_car ) * GlobalC::ucell.lat0;
+                            ModuleBase::Vector3<double> r_distance = 
+                                ( GlobalC::ucell.atoms[it2].tau[ia2] - 
+                                  GlobalC::ucell.atoms[it1].tau[ia1] + R_car ) 
+                                * GlobalC::ucell.lat0;
+
+                            double overlap_o = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( 
+                                   origin_point, r_distance, im1, im2 );
+
+                            double overlap_x = -1 * factor * 
+                                   center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( 
+                                   origin_point, r_distance, im1, 1, im2 ); // m =  1
 
-                            double overlap_o = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( origin_point, r_distance, im1, im2 );
-                            double overlap_x = -1 * factor * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( origin_point, r_distance, im1, 1, im2 ); // m =  1
-                            double overlap_y = -1 * factor * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( origin_point, r_distance, im1, 2, im2 ); // m = -1
-                            double overlap_z =      factor * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( origin_point, r_distance, im1, 0, im2 ); // m =  0
-                            ModuleBase::Vector3<double> temp_prp = ModuleBase::Vector3<double>(overlap_x, overlap_y, overlap_z) + GlobalC::ucell.atoms[it1].tau[ia1] * GlobalC::ucell.lat0 * overlap_o;
+                            double overlap_y = -1 * factor * 
+                                   center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( 
+                                   origin_point, r_distance, im1, 2, im2 ); // m = -1
+
+                            double overlap_z = factor * 
+                                   center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap( 
+                                   origin_point, r_distance, im1, 0, im2 ); // m =  0
+
+                            ModuleBase::Vector3<double> temp_prp = 
+                            ModuleBase::Vector3<double>(overlap_x, overlap_y, overlap_z) 
+                            + GlobalC::ucell.atoms[it1].tau[ia1] * GlobalC::ucell.lat0 * overlap_o;
 
                             if (std::abs(temp_prp.x) > sparse_threshold)
                             {
@@ -642,4 +683,4 @@ void cal_r_overlap_R::out_rR_other(const int &istep, const std::set<Abfs::Vector
 
     ModuleBase::timer::tick("cal_r_overlap_R", "out_rR_other");
     return;
-}
\ No newline at end of file
+}
diff --git a/source/module_io/output_mat_sparse.cpp b/source/module_io/output_mat_sparse.cpp
index 0c82f71a350..d68e218eb2e 100644
--- a/source/module_io/output_mat_sparse.cpp
+++ b/source/module_io/output_mat_sparse.cpp
@@ -1,5 +1,4 @@
 #include "output_mat_sparse.h"
-
 #include "cal_r_overlap_R.h"
 #include "write_HS_R.h"
 
@@ -16,6 +15,7 @@ namespace ModuleIO
         const ModuleBase::matrix& v_eff,
         const Parallel_Orbitals& pv,
         LCAO_Hamilt& UHM,
+        LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
         Gint_k& gint_k, // mohan add 2024-04-01
         LCAO_Matrix& LM,
         const K_Vectors& kv,
@@ -28,6 +28,7 @@ namespace ModuleIO
       _v_eff(v_eff),
       _pv(pv),
       _UHM(UHM),
+      _gen_h(gen_h), // mohan add 2024-04-02
       _gint_k(gint_k), // mohan add 2024-04-01
       _LM(LM),
       _kv(kv),
@@ -44,22 +45,26 @@ void Output_Mat_Sparse<double>::write(void)
 template<>
 void Output_Mat_Sparse<std::complex<double>>::write(void)
 {
+    //! generate a file containing the Hamiltonian and S(overlap) matrices 
     if (_out_mat_hsR)
     {
         output_HS_R(_istep, this->_v_eff, this->_UHM, _kv, _p_ham);
     }
 
+    //! generate a file containing the kinetic energy matrix
     if (_out_mat_t)
     {
-        output_T_R(_istep, this->_UHM); // LiuXh add 2019-07-15
+        output_T_R(_istep, this->_UHM, this->_gen_h); // LiuXh add 2019-07-15
     }
 
+    //! generate a file containing the derivatives of the Hamiltonian matrix (in Ry/Bohr)
     if (_out_mat_dh)
     {
 		output_dH_R(
 				_istep, 
 				this->_v_eff, 
 				this->_UHM, 
+                this->_gen_h,
 				this->_gint_k, // mohan add 2024-04-01
 				this->_LM,
 				_kv); // LiuXh add 2019-07-15
diff --git a/source/module_io/output_mat_sparse.h b/source/module_io/output_mat_sparse.h
index 8012b108bca..c655226089b 100644
--- a/source/module_io/output_mat_sparse.h
+++ b/source/module_io/output_mat_sparse.h
@@ -23,6 +23,7 @@ namespace ModuleIO
         const ModuleBase::matrix& v_eff,
         const Parallel_Orbitals& pv,
         LCAO_Hamilt& UHM,
+        LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
         Gint_k &gint_k, // mohan add 2024-04-01
         LCAO_Matrix& LM,
         const K_Vectors& kv,
@@ -31,19 +32,38 @@ namespace ModuleIO
     void write() override;
 
   private:
+
+    //! generate a file containing the Hamiltonian and S(overlap) matrices 
     int _out_mat_hsR;
+
+    //! generate a file containing the derivatives of the Hamiltonian matrix (in Ry/Bohr)
     int _out_mat_dh;
+
+    //! generate a file containing the kinetic energy matrix
     int _out_mat_t;
+
+    //! generate a file containing the matrix representation of the position matrix (in Bohr)
     int _out_mat_r;
+
     int _istep;
+
     const ModuleBase::matrix& _v_eff;
+
     const Parallel_Orbitals& _pv;
+
     LCAO_Hamilt& _UHM;
+
+    LCAO_gen_fixedH& _gen_h; // mohan add 2024-04-02
+
     Gint_k& _gint_k; // mohan add 2024-04-01
+
     LCAO_Matrix& _LM;
+
     const K_Vectors& _kv;
+
     hamilt::Hamilt<T>* _p_ham;
 };
+
 } // namespace ModuleIO
 
 #endif // OUTPUT_MAT_SPARSE_H
diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp
index 182fdba02a0..d6cb45fb54b 100644
--- a/source/module_io/td_current_io.cpp
+++ b/source/module_io/td_current_io.cpp
@@ -11,40 +11,48 @@
 
 #ifdef __LCAO
 //init DSloc_R for current calculation
-void ModuleIO::Init_DS_tmp(const Parallel_Orbitals& pv,
-                              LCAO_Hamilt& UHM)
+void ModuleIO::Init_DS_tmp(
+		const Parallel_Orbitals& pv,
+		LCAO_Matrix &lm,
+		LCAO_gen_fixedH &gen_h)
 {    
     ModuleBase::TITLE("ModuleIO", "Init_DS_tmp");
     ModuleBase::timer::tick("ModuleIO", "Init_DS_tmp");
     const int nnr = pv.nnr;
-    UHM.LM->DSloc_Rx = new double[nnr];
-    UHM.LM->DSloc_Ry = new double[nnr];
-    UHM.LM->DSloc_Rz = new double[nnr];
-    const auto init_DSloc_Rxyz = [&UHM, nnr](int num_threads, int thread_id) {
+    lm.DSloc_Rx = new double[nnr];
+    lm.DSloc_Ry = new double[nnr];
+    lm.DSloc_Rz = new double[nnr];
+
+    const auto init_DSloc_Rxyz = [&lm, nnr](int num_threads, int thread_id) {
         int beg, len;
         ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, nnr, 1024, beg, len);
-        ModuleBase::GlobalFunc::ZEROS(UHM.LM->DSloc_Rx + beg, len);
-        ModuleBase::GlobalFunc::ZEROS(UHM.LM->DSloc_Ry + beg, len);
-        ModuleBase::GlobalFunc::ZEROS(UHM.LM->DSloc_Rz + beg, len);
+        ModuleBase::GlobalFunc::ZEROS(lm.DSloc_Rx + beg, len);
+        ModuleBase::GlobalFunc::ZEROS(lm.DSloc_Ry + beg, len);
+        ModuleBase::GlobalFunc::ZEROS(lm.DSloc_Rz + beg, len);
     };
+
     ModuleBase::OMP_PARALLEL(init_DSloc_Rxyz);
     bool cal_deri = true;
-    UHM.genH.build_ST_new('S', cal_deri, GlobalC::ucell, UHM.genH.LM->SlocR.data());
+    gen_h.build_ST_new('S', cal_deri, GlobalC::ucell, lm.SlocR.data());
 
     ModuleBase::timer::tick("ModuleIO", "Init_DS_tmp");
     return;
 }
+
+
 //destory DSloc_R so it can be used normally in the following force calculation
-void ModuleIO::destory_DS_tmp(LCAO_Hamilt& UHM)
+void ModuleIO::destory_DS_tmp(LCAO_Matrix &lm)
 {
     ModuleBase::TITLE("ModuleIO", "destory_DS_tmp");
     ModuleBase::timer::tick("ModuleIO", "destory_DS_tmp");
-    delete[] UHM.LM->DSloc_Rx;
-    delete[] UHM.LM->DSloc_Ry;
-    delete[] UHM.LM->DSloc_Rz;
+    delete[] lm.DSloc_Rx;
+    delete[] lm.DSloc_Ry;
+    delete[] lm.DSloc_Rz;
     ModuleBase::timer::tick("ModuleIO", "destory_DS_tmp");
     return;
 }
+
+
 void ModuleIO::cal_tmp_DM(elecstate::DensityMatrix<std::complex<double>, double>& DM, const int ik, const int nspin)
 {
     ModuleBase::TITLE("ModuleIO", "cal_tmp_DM");
@@ -133,14 +141,15 @@ void ModuleIO::write_current(const int istep,
                                 const elecstate::ElecState* pelec,
                                 const K_Vectors& kv,
                                 const Parallel_Orbitals* pv,
-                                Record_adj& ra,
-                                LCAO_Hamilt& UHM)
+								Record_adj& ra,
+								LCAO_Matrix &lm, // mohan add 2024-04-02
+								LCAO_gen_fixedH &gen_h) // mohan add 2024-04-02
 {
 
     ModuleBase::TITLE("ModuleIO", "write_current");
     ModuleBase::timer::tick("ModuleIO", "write_current");
     //Init_DS_tmp
-    Init_DS_tmp(*pv, UHM);
+    Init_DS_tmp(*pv, lm, gen_h);
     // construct a DensityMatrix object
     elecstate::DensityMatrix<std::complex<double>, double> DM(&kv,pv,GlobalV::NSPIN);
     
@@ -209,9 +218,9 @@ void ModuleIO::write_current(const int istep,
                             // here do not sum over spin due to EDM.sum_DMR_spin();
                             double edm2d1 = tmp_matrix->get_value(mu,nu);
                             double edm2d2 = 2.0 * edm2d1;
-                            current_ik[0] -= edm2d2 * UHM.LM->DSloc_Rx[irr];
-                            current_ik[1] -= edm2d2 * UHM.LM->DSloc_Ry[irr];
-                            current_ik[2] -= edm2d2 * UHM.LM->DSloc_Rz[irr];
+                            current_ik[0] -= edm2d2 * lm.DSloc_Rx[irr];
+                            current_ik[1] -= edm2d2 * lm.DSloc_Ry[irr];
+                            current_ik[2] -= edm2d2 * lm.DSloc_Rz[irr];
                             ++local_total_irr;
                             ++irr;
                         } // end kk
@@ -241,7 +250,7 @@ void ModuleIO::write_current(const int istep,
         //write end
         ModuleBase::timer::tick("ModuleIO", "write_current");
     }//end nks
-    destory_DS_tmp(UHM);
+    destory_DS_tmp(lm);
     return;
 }
 #endif //__LCAO
diff --git a/source/module_io/td_current_io.h b/source/module_io/td_current_io.h
index 65a9080b096..63643a179e7 100644
--- a/source/module_io/td_current_io.h
+++ b/source/module_io/td_current_io.h
@@ -15,17 +15,21 @@ void write_current(const int istep,
                     const elecstate::ElecState* pelec,
                     const K_Vectors& kv,
                     const Parallel_Orbitals* pv,
-                    Record_adj& ra,
-                    LCAO_Hamilt& UHM);
+					Record_adj& ra,
+					LCAO_Matrix &lm, // mohan add 2024-04-02
+					LCAO_gen_fixedH &gen_h); // mohan add 2024-04-02
 
 /// @brief calculate sum_n[𝜌_(𝑛𝑘,𝜇𝜈)] for current calculation
 void cal_tmp_DM(elecstate::DensityMatrix<std::complex<double>, double>& DM, const int ik, const int nspin);
 
 /// @brief Init DS_locR for currrent calculation
-void Init_DS_tmp(const Parallel_Orbitals& pv,LCAO_Hamilt& UHM);
+void Init_DS_tmp(
+		const Parallel_Orbitals& pv,
+		LCAO_Matrix &lm,
+		LCAO_gen_fixedH &gen_h);
 
 /// @brief DS_locR will be initialized again in force calculation, so it must be destoryed here.
-void destory_DS_tmp(LCAO_Hamilt& UHM);
+void destory_DS_tmp(LCAO_Matrix &lm);
 
 #endif // __LCAO
 }
diff --git a/source/module_io/write_HS_R.cpp b/source/module_io/write_HS_R.cpp
index 5869b6ee5f3..41b19f1d0f1 100644
--- a/source/module_io/write_HS_R.cpp
+++ b/source/module_io/write_HS_R.cpp
@@ -3,12 +3,15 @@
 #include "module_base/timer.h"
 #include "write_HS_sparse.h"
 
+#include "module_hamilt_lcao/hamilt_lcaodft/sparse_format.h"
+
+
 // if 'binary=true', output binary file.
 // The 'sparse_threshold' is the accuracy of the sparse matrix. 
 // If the absolute value of the matrix element is less than or equal to the 'sparse_threshold', it will be ignored.
 void ModuleIO::output_HS_R(const int& istep,
                            const ModuleBase::matrix& v_eff,
-                           LCAO_Hamilt& UHM,
+                           LCAO_Hamilt& uhm,
                            const K_Vectors& kv,
                            hamilt::Hamilt<std::complex<double>>* p_ham,
                            const std::string& SR_filename,
@@ -22,13 +25,15 @@ void ModuleIO::output_HS_R(const int& istep,
 
     if(GlobalV::NSPIN==1||GlobalV::NSPIN==4)
     {
+        // mohan add 2024-04-02
+        const int spin_now = 0;
         // jingan add 2021-6-4, modify 2021-12-2
-        UHM.cal_HSR_sparse(0, sparse_threshold, kv.nmp, p_ham);
+        uhm.cal_HSR_sparse(spin_now, sparse_threshold, kv.nmp, p_ham);
     }
     else if(GlobalV::NSPIN==2)
     {
         // save HR of current_spin first
-        UHM.cal_HSR_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, kv.nmp, p_ham);
+        uhm.cal_HSR_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, kv.nmp, p_ham);
         // cal HR of the other spin
         if(GlobalV::VL_IN_H)
         {
@@ -46,11 +51,11 @@ void ModuleIO::output_HS_R(const int& istep,
             p_ham->refresh();
             p_ham->updateHk(ik);
         }
-        UHM.cal_HSR_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, kv.nmp, p_ham);
+        uhm.cal_HSR_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, kv.nmp, p_ham);
     }
 
-    ModuleIO::save_HSR_sparse(istep, *UHM.LM, sparse_threshold, binary, SR_filename, HR_filename_up, HR_filename_down);
-    UHM.destroy_all_HSR_sparse();
+    ModuleIO::save_HSR_sparse(istep, *uhm.LM, sparse_threshold, binary, SR_filename, HR_filename_up, HR_filename_down);
+    uhm.destroy_all_HSR_sparse();
 
     ModuleBase::timer::tick("ModuleIO","output_HS_R"); 
     return;
@@ -58,7 +63,8 @@ void ModuleIO::output_HS_R(const int& istep,
 
 void ModuleIO::output_dH_R(const int& istep,
                            const ModuleBase::matrix& v_eff,
-                           LCAO_Hamilt& UHM,
+                           LCAO_Hamilt& uhm,
+                           LCAO_gen_fixedH& gen_h, // mohan add 2024-04-02
                            Gint_k& gint_k,  // mohan add 2024-04-01
                            LCAO_Matrix &lm,  // mohan add 2024-04-01
                            const K_Vectors& kv,
@@ -74,8 +80,14 @@ void ModuleIO::output_dH_R(const int& istep,
     {
         // mohan add 2024-04-01
         assert(GlobalV::CURRENT_SPIN==0);
-        UHM.cal_dH_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, gint_k);
-    }
+
+		sparse_format::cal_dH(
+                lm,
+				gen_h,
+				GlobalV::CURRENT_SPIN, 
+				sparse_threshold, 
+				gint_k);
+	}
     else if(GlobalV::NSPIN==2)
     {
         for (int ik = 0; ik < kv.nks; ik++)
@@ -99,14 +111,19 @@ void ModuleIO::output_dH_R(const int& istep,
                     }
                 }
 
-                UHM.cal_dH_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, gint_k);
-            }
+				sparse_format::cal_dH(
+                        lm,
+						gen_h,
+						GlobalV::CURRENT_SPIN, 
+						sparse_threshold, 
+						gint_k);
+			}
         }
     }
 
     // mohan update 2024-04-01
     ModuleIO::save_dH_sparse(istep, lm, sparse_threshold, binary);
-    UHM.destroy_dH_R_sparse();
+    uhm.destroy_dH_R_sparse();
 
     gint_k.destroy_pvdpR();
 
@@ -115,7 +132,7 @@ void ModuleIO::output_dH_R(const int& istep,
 }
 
 void ModuleIO::output_S_R(
-    LCAO_Hamilt &UHM,
+    LCAO_Hamilt &uhm,
     hamilt::Hamilt<std::complex<double>>* p_ham,
     const std::string &SR_filename,
     const bool &binary,
@@ -124,9 +141,20 @@ void ModuleIO::output_S_R(
     ModuleBase::TITLE("ModuleIO","output_S_R");
     ModuleBase::timer::tick("ModuleIO","output_S_R"); 
 
-    UHM.cal_SR_sparse(sparse_threshold, p_ham);
-    ModuleIO::save_sparse(UHM.LM->SR_sparse, UHM.LM->all_R_coor, sparse_threshold, binary, SR_filename, *UHM.LM->ParaV, "S", 0);
-    UHM.destroy_all_HSR_sparse();
+    uhm.cal_SR_sparse(sparse_threshold, p_ham);
+
+	ModuleIO::save_sparse(
+			uhm.LM->SR_sparse, 
+			uhm.LM->all_R_coor,
+			sparse_threshold, 
+			binary, 
+			SR_filename,
+			*uhm.LM->ParaV, 
+			"S", 
+			0
+			);
+
+    uhm.destroy_all_HSR_sparse();
 
     ModuleBase::timer::tick("ModuleIO","output_S_R");
     return;
@@ -134,7 +162,8 @@ void ModuleIO::output_S_R(
 
 void ModuleIO::output_T_R(
     const int istep,
-    LCAO_Hamilt &UHM,
+    LCAO_Hamilt &uhm,
+    LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
     const std::string &TR_filename,
     const bool &binary,
     const double &sparse_threshold
@@ -153,9 +182,20 @@ void ModuleIO::output_T_R(
         sst << GlobalV::global_out_dir << TR_filename;
     }
 
-    UHM.cal_TR_sparse(sparse_threshold);
-    ModuleIO::save_sparse(UHM.LM->TR_sparse, UHM.LM->all_R_coor, sparse_threshold, binary, sst.str().c_str(), *UHM.LM->ParaV, "T", istep);
-    UHM.destroy_TR_sparse();
+    uhm.cal_TR_sparse(gen_h, sparse_threshold);
+
+	ModuleIO::save_sparse(
+			uhm.LM->TR_sparse, 
+			uhm.LM->all_R_coor,
+			sparse_threshold, 
+			binary, 
+		    sst.str().c_str(),
+			*uhm.LM->ParaV, 
+			"T", 
+			istep
+			);
+
+    uhm.destroy_TR_sparse();
 
     ModuleBase::timer::tick("ModuleIO","output_T_R");
     return;
diff --git a/source/module_io/write_HS_R.h b/source/module_io/write_HS_R.h
index 738af1d3a7d..b61b3daf7f3 100644
--- a/source/module_io/write_HS_R.h
+++ b/source/module_io/write_HS_R.h
@@ -3,8 +3,8 @@
 
 #include "module_base/matrix.h"
 #include "module_cell/klist.h"
-#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h"
 #include "module_hamilt_general/hamilt.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h"
 
 namespace ModuleIO
 {
@@ -22,8 +22,9 @@ namespace ModuleIO
 
         void output_dH_R(
             const int &istep,
-            const ModuleBase::matrix& v_eff,
-			LCAO_Hamilt &UHM,
+			const ModuleBase::matrix& v_eff,
+			LCAO_Hamilt &uhm,
+			LCAO_gen_fixedH& gen_h, // mohan add 2024-04-02
 			Gint_k& gint_k,  // mohan add 2024-04-01
 			LCAO_Matrix &lm,  // mohan add 2024-04-01
             const K_Vectors& kv,
@@ -33,6 +34,7 @@ namespace ModuleIO
         void output_T_R(
             const int istep,
             LCAO_Hamilt &UHM,
+            LCAO_gen_fixedH &gen_h, // mohan add 2024-04-02
             const std::string& TR_filename = "data-TR-sparse_SPIN0.csr",
             const bool& binary = false,
             const double& sparse_threshold = 1e-10);
diff --git a/source/module_io/write_HS_sparse.cpp b/source/module_io/write_HS_sparse.cpp
index 4abd1b6acf5..bfd534585a6 100644
--- a/source/module_io/write_HS_sparse.cpp
+++ b/source/module_io/write_HS_sparse.cpp
@@ -338,6 +338,7 @@ void ModuleIO::save_HSR_sparse(
     return;
 }
 
+
 void ModuleIO::save_dH_sparse(
     const int &istep,
     LCAO_Matrix &lm,
@@ -598,7 +599,7 @@ void ModuleIO::save_dH_sparse(
                     g1z[ispin].write(reinterpret_cast<char *>(&dRx), sizeof(int));
                     g1z[ispin].write(reinterpret_cast<char *>(&dRy), sizeof(int));
                     g1z[ispin].write(reinterpret_cast<char *>(&dRz), sizeof(int));
-                    g1z[ispin].write(reinterpret_cast<char *>(&dHz_nonzero_num[ispin][count]), sizeof(int));                                        
+                    g1z[ispin].write(reinterpret_cast<char *>(&dHz_nonzero_num[ispin][count]), sizeof(int)); 
                 }
             }
             else
@@ -655,9 +656,18 @@ void ModuleIO::save_dH_sparse(
 
     if(GlobalV::DRANK==0) 
     {
-        for (int ispin = 0; ispin < spin_loop; ++ispin) g1x[ispin].close();
-        for (int ispin = 0; ispin < spin_loop; ++ispin) g1y[ispin].close();
-        for (int ispin = 0; ispin < spin_loop; ++ispin) g1z[ispin].close();
+		for (int ispin = 0; ispin < spin_loop; ++ispin) 
+		{
+			g1x[ispin].close();
+		}
+		for (int ispin = 0; ispin < spin_loop; ++ispin) 
+		{
+			g1y[ispin].close();
+		}
+		for (int ispin = 0; ispin < spin_loop; ++ispin) 
+		{
+			g1z[ispin].close();
+		}
     }
     
     for (int ispin = 0; ispin < spin_loop; ++ispin) 
@@ -695,16 +705,28 @@ void ModuleIO::save_sparse(
     for (auto& R_coor : all_R_coor)
     {
         auto iter = smat.find(R_coor);
-        if (iter != smat.end())
-            for (auto& row_loop : iter->second)
-                nonzero_num[count] += row_loop.second.size();
+		if (iter != smat.end())
+		{
+			for (auto& row_loop : iter->second)
+			{
+				nonzero_num[count] += row_loop.second.size();
+			}
+		}
         ++count;
     }
-    if (reduce)Parallel_Reduce::reduce_all(nonzero_num.data(), total_R_num);
+	if (reduce)
+	{
+		Parallel_Reduce::reduce_all(nonzero_num.data(), total_R_num);
+	}
 
     int output_R_number = 0;
-    for (int index = 0; index < total_R_num; ++index)
-        if (nonzero_num[index] != 0) ++output_R_number;
+	for (int index = 0; index < total_R_num; ++index)
+	{
+		if (nonzero_num[index] != 0) 
+		{
+			++output_R_number;
+		}
+	}
 
     std::stringstream sss;
     sss << filename;
@@ -713,10 +735,14 @@ void ModuleIO::save_sparse(
     {
         if (binary)
         {
-            if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && istep)
-                ofs.open(sss.str().c_str(), std::ios::binary | std::ios::app);
-            else
-                ofs.open(sss.str().c_str(), std::ios::binary);
+			if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && istep)
+			{
+				ofs.open(sss.str().c_str(), std::ios::binary | std::ios::app);
+			}
+			else
+			{
+				ofs.open(sss.str().c_str(), std::ios::binary);
+			}
             ofs.write(reinterpret_cast<char*>(0), sizeof(int));
             ofs.write(reinterpret_cast<char*>(&GlobalV::NLOCAL), sizeof(int));
             ofs.write(reinterpret_cast<char*>(&output_R_number), sizeof(int));
@@ -724,9 +750,13 @@ void ModuleIO::save_sparse(
         else
         {
             if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && istep)
-                ofs.open(sss.str().c_str(), std::ios::app);
-            else
-                ofs.open(sss.str().c_str());
+			{
+				ofs.open(sss.str().c_str(), std::ios::app);
+			}
+			else
+			{
+				ofs.open(sss.str().c_str());
+			}
             ofs << "STEP: " << std::max(istep, 0) << std::endl;
             ofs << "Matrix Dimension of " + label + "(R): " << GlobalV::NLOCAL << std::endl;
             ofs << "Matrix number of " + label + "(R): " << output_R_number << std::endl;
@@ -764,7 +794,10 @@ void ModuleIO::save_sparse(
         output_single_R(ofs, smat.at(R_coor), sparse_threshold, binary, pv, reduce);
         ++count;
     }
-    if (!reduce || GlobalV::DRANK == 0) ofs.close();
+	if (!reduce || GlobalV::DRANK == 0) 
+	{
+		ofs.close();
+	}
 
     ModuleBase::timer::tick("ModuleIO", "save_sparse");
 }
@@ -779,6 +812,8 @@ template void ModuleIO::save_sparse<double>(
     const std::string&,
     const int&,
     const bool&);
+
+
 template void ModuleIO::save_sparse<std::complex<double>>(
     const std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>>&,
     const std::set<Abfs::Vector3_Order<int>>&,
@@ -788,4 +823,4 @@ template void ModuleIO::save_sparse<std::complex<double>>(
     const Parallel_Orbitals&,
     const std::string&,
     const int&,
-    const bool&);
\ No newline at end of file
+    const bool&);