diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
index 419dc48ca3e..de479128f3a 100644
--- a/source/module_esolver/esolver_ks.cpp
+++ b/source/module_esolver/esolver_ks.cpp
@@ -392,7 +392,9 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
 		ModuleBase::timer::tick(this->classname, "run");
 
 		this->before_scf(istep); //Something else to do before the iter loop
+
 		ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
+
 		if(this->maxniter > 0)  
 		{
 			this->print_head(); //print the headline on the screen.
@@ -499,9 +501,6 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
                (std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::system_clock::now() 
                 - iterstart)).count() / static_cast<double>(1e6);
 #endif
-			/*
-			   SCF print: G1    -3.435545e+03  0.000000e+00   3.607e-01  2.862e-01
-			 */
 
 			double dkin = 0.0; // for meta-GGA
 			if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
@@ -525,7 +524,10 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
 			{
 				this->niter = iter;
 				bool stop = this->do_after_converge(iter);
-				if(stop) break;
+				if(stop) 
+				{
+					break;
+				}
 			}
 
 			// notice for restart
@@ -544,7 +546,9 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
 				this->conv_elec
 				);
 #endif //__RAPIDJSON 
-		after_scf(istep);
+
+		this->after_scf(istep);
+
 		ModuleBase::timer::tick(this->classname, "run");
 	} 
 
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
index 03433f9198b..3379777f905 100644
--- a/source/module_esolver/esolver_ks_lcao.cpp
+++ b/source/module_esolver/esolver_ks_lcao.cpp
@@ -1,4 +1,5 @@
 #include "esolver_ks_lcao.h"
+#include "module_base/tool_title.h"
 
 #include "module_base/global_variable.h"
 #include "module_io/dos_nao.h"
@@ -21,7 +22,6 @@
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/print_info.h"
 #ifdef __EXX
-// #include "module_rpa/rpa.h"
 #include "module_ri/RPA_LRI.h"
 #endif
 
@@ -81,6 +81,8 @@ template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
 {
     ModuleBase::TITLE("ESolver_KS_LCAO", "init");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "init");
+
     // if we are only calculating S, then there is no need
     // to prepare for potentials and so on
 
@@ -123,7 +125,7 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
     // Init Basis should be put outside of Ensolver.
     // * reading the localized orbitals/projectors
     // * construct the interpolation tables.
-    this->Init_Basis_lcao(this->orb_con, inp, ucell);
+    this->init_basis_lcao(this->orb_con, inp, ucell);
     //------------------init Basis_lcao----------------------
 
     //! pass Hamilt-pointer to Operator
@@ -137,6 +139,7 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
 
     if (GlobalV::CALCULATION == "get_S")
     {
+        ModuleBase::timer::tick("ESolver_KS_LCAO", "init");
         return;
     }
 
@@ -230,13 +233,19 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
     if (GlobalV::ocp)
     {
         this->pelec->fixed_weights(GlobalV::ocp_kb);
-    }
+	}
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "init");
+	return;
 }
 
 
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::init_after_vc(Input& inp, UnitCell& ucell)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "init_after_vc");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc");
+
 	ESolver_KS<TK>::init_after_vc(inp, ucell);
 
     if (GlobalV::md_prec_level == 2)
@@ -269,13 +278,18 @@ void ESolver_KS_LCAO<TK, TR>::init_after_vc(Input& inp, UnitCell& ucell)
                                                         &(this->pelec->f_en.etxc),
                                                         &(this->pelec->f_en.vtxc));
         }
-    }
+	}
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc");
+	return;
 }
 
 
 template <typename TK, typename TR>
 double ESolver_KS_LCAO<TK, TR>::cal_energy()
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "cal_energy");
+
     return this->pelec->f_en.etot;
 }
 
@@ -283,12 +297,19 @@ double ESolver_KS_LCAO<TK, TR>::cal_energy()
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::cal_force(ModuleBase::matrix& force)
 {
-	Force_Stress_LCAO<TK> FSL(this->RA, GlobalC::ucell.nat);
-	FSL.getForceStress(GlobalV::CAL_FORCE,
+    ModuleBase::TITLE("ESolver_KS_LCAO", "cal_force");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_force");
+
+	Force_Stress_LCAO<TK> fsl(this->RA, GlobalC::ucell.nat);
+
+	fsl.getForceStress(
+            GlobalV::CAL_FORCE,
 			GlobalV::CAL_STRESS,
 			GlobalV::TEST_FORCE,
 			GlobalV::TEST_STRESS,
 			this->LOC,
+            this->orb_con.ParaV, 
+            this->LM,
 			this->pelec,
 			this->psi,
 			this->UHM,
@@ -308,12 +329,17 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(ModuleBase::matrix& force)
 	this->RA.delete_grid();
 
 	this->have_force = true;
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_force");
 }
 
 
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::cal_stress(ModuleBase::matrix& stress)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "cal_stress");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress");
+
     if (!this->have_force)
     {
         ModuleBase::matrix fcs;
@@ -321,12 +347,17 @@ void ESolver_KS_LCAO<TK, TR>::cal_stress(ModuleBase::matrix& stress)
     }
     stress = this->scs; // copy the stress
     this->have_force = false;
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress");
 }
 
 
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::post_process(void)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "post_process");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "post_process");
+
     GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl;
     GlobalV::ofs_running << std::setprecision(16);
     GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
@@ -413,15 +444,18 @@ void ESolver_KS_LCAO<TK, TR>::post_process(void)
 				GlobalV::NBANDS,
 				this->p_hamilt);
 	}
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "post_process");
 }
 
 
 template <typename TK, typename TR>
-void ESolver_KS_LCAO<TK, TR>::Init_Basis_lcao(
+void ESolver_KS_LCAO<TK, TR>::init_basis_lcao(
 		ORB_control& orb_con, 
 		Input& inp, 
 		UnitCell& ucell)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "init_basis_lcao");
+
     // autoset NB2D first
     if (GlobalV::NB2D == 0)
     {
@@ -500,12 +534,16 @@ void ESolver_KS_LCAO<TK, TR>::Init_Basis_lcao(
         this->orb_con.setup_2d_division(GlobalV::ofs_running, GlobalV::ofs_warning);
         this->orb_con.ParaV.set_atomic_trace(GlobalC::ucell.get_iat2iwt(), GlobalC::ucell.nat, GlobalV::NLOCAL);
     }
+
+    return;
 }
 
 
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "iter_init");
+
     if (iter == 1)
     {
         this->p_chgmix->init_mixing(); // init mixing
@@ -645,6 +683,8 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density");
+
     // save input rho
     this->pelec->charge->save_rho_before_sum_band();
     // save density matrix for mixing
@@ -756,6 +796,8 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "update_pot");
+
     // print Hamiltonian and Overlap matrix
     if (this->conv_elec)
     {
@@ -852,6 +894,8 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::iter_finish(int iter)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish");
+
     // mix density matrix
     if (GlobalV::MIXING_RESTART > 0 && iter >= this->p_chgmix->mixing_restart && GlobalV::MIXING_DMR )
     {
@@ -930,6 +974,8 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int iter)
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf");
+
     // save charge difference into files for charge extrapolation
     if (GlobalV::CALCULATION != "scf")
     {
@@ -1015,6 +1061,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
             this->LOWF.ParaV,
             *(this->psi),
             dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM());
+
     ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels");
 
 #endif
@@ -1042,6 +1089,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
         } // qifeng add 2019/9/10, jiyy modify 2023/2/27, liuyu move here 2023-04-18
     }
 
+    // spin constrain calculations, added by Tianqi Zhao.
     if (GlobalV::sc_mag_switch)
     {
         SpinConstrain<TK, psi::DEVICE_CPU>& sc = SpinConstrain<TK, psi::DEVICE_CPU>::getScInstance();
@@ -1052,6 +1100,8 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
     {
         RA.delete_grid();
     }
+
+    // quasi-orbitals, added by Yike Huang.
     if(GlobalV::qo_switch)
     {
         toQO tqo(GlobalV::qo_basis, GlobalV::qo_strategy, GlobalV::qo_thr, GlobalV::qo_screening_coeff);
@@ -1069,6 +1119,8 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
 template <typename TK, typename TR>
 bool ESolver_KS_LCAO<TK, TR>::do_after_converge(int& iter)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO","do_after_converge");
+
 #ifdef __EXX
     if (GlobalC::exx_info.info_ri.real_number)
     {
@@ -1088,11 +1140,13 @@ bool ESolver_KS_LCAO<TK, TR>::do_after_converge(int& iter)
 				iter);
 	}
 #endif // __EXX
+
     if(GlobalV::dft_plus_u)
     {
         // use the converged occupation matrix for next MD/Relax SCF calculation
         GlobalC::dftu.initialed_locale = true;
     }
+
     return true;
 }
 
@@ -1127,17 +1181,17 @@ ModuleIO::Output_DM1 ESolver_KS_LCAO<TK, TR>::create_Output_DM1(int istep)
 template <typename TK, typename TR>
 ModuleIO::Output_Mat_Sparse<TK> ESolver_KS_LCAO<TK, TR>::create_Output_Mat_Sparse(int istep)
 {
-        return ModuleIO::Output_Mat_Sparse<TK>(hsolver::HSolverLCAO<TK>::out_mat_hsR,
-            hsolver::HSolverLCAO<TK>::out_mat_dh,
-            hsolver::HSolverLCAO<TK>::out_mat_t,
-            INPUT.out_mat_r,
-            istep,
-            this->pelec->pot->get_effective_v(),
-            *this->LOWF.ParaV,
-            this->UHM,
-            this->LM,
-            this->kv,
-            this->p_hamilt);
+	return ModuleIO::Output_Mat_Sparse<TK>(hsolver::HSolverLCAO<TK>::out_mat_hsR,
+			hsolver::HSolverLCAO<TK>::out_mat_dh,
+			hsolver::HSolverLCAO<TK>::out_mat_t,
+			INPUT.out_mat_r,
+			istep,
+			this->pelec->pot->get_effective_v(),
+			*this->LOWF.ParaV,
+			this->UHM,
+			this->LM,
+			this->kv,
+			this->p_hamilt);
 }
 
 
diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h
index 269bb2b3a1e..49d95cdfe55 100644
--- a/source/module_esolver/esolver_ks_lcao.h
+++ b/source/module_esolver/esolver_ks_lcao.h
@@ -32,9 +32,9 @@ namespace ModuleESolver
 
         double cal_energy() override;
 
-        void cal_force(ModuleBase::matrix& force) override;
+        void cal_force(ModuleBase::matrix &force) override;
 
-        void cal_stress(ModuleBase::matrix& stress) override;
+        void cal_stress(ModuleBase::matrix &stress) override;
 
         void post_process() override;
 
@@ -89,7 +89,7 @@ namespace ModuleESolver
         ModuleBase::matrix scs;
         bool have_force = false;
 
-        void Init_Basis_lcao(ORB_control& orb_con, Input& inp, UnitCell& ucell);
+        void init_basis_lcao(ORB_control& orb_con, Input& inp, UnitCell& ucell);
 
         //--------------common for all calculation, not only scf-------------
         // set matrix and grid integral
diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp
index a4ff74d1501..a9bd2e44313 100644
--- a/source/module_esolver/esolver_ks_lcao_tddft.cpp
+++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp
@@ -90,7 +90,7 @@ void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell)
     // Init Basis should be put outside of Ensolver.
     // * reading the localized orbitals/projectors
     // * construct the interpolation tables.
-    this->Init_Basis_lcao(this->orb_con, inp, ucell);
+    this->init_basis_lcao(this->orb_con, inp, ucell);
     //------------------init Basis_lcao----------------------
 
     //------------------init Hamilt_lcao----------------------
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
index 3608202ee8d..4fa805e78de 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
@@ -29,7 +29,9 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                                           const bool isstress,
                                           const bool istestf,
                                           const bool istests,
-                                          Local_Orbital_Charge& loc,
+										  Local_Orbital_Charge &loc,
+										  Parallel_Orbitals &pv,
+										  LCAO_Matrix &lm,
                                           const elecstate::ElecState* pelec,
                                           const psi::Psi<T>* psi,
                                           LCAO_Hamilt& uhm,
@@ -155,7 +157,9 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
 #else
                                      svl_dphi,
 #endif
-                                     uhm,
+									 uhm,
+									 pv,
+									 lm,
                                      kv);
     // implement vdw force or stress here
     //  Peize Lin add 2014-04-04, update 2021-03-09
@@ -242,10 +246,12 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
             tmp_dftu.cal_force_stress(isforce, isstress, force_dftu, stress_dftu);
         }
     }
+
     if (!GlobalV::GAMMA_ONLY_LOCAL)
     {
-        this->flk.finish_k();
+        this->flk.finish_k(lm);
     }
+
 #ifdef __EXX
     // Force and Stress contribution from exx
     ModuleBase::matrix force_exx;
@@ -712,7 +718,8 @@ void Force_Stress_LCAO<T>::calForcePwPart(ModuleBase::matrix& fvl_dvl,
 
 // overlap, kinetic, nonlocal pseudopotential, Local potential terms in force and stress
 template<>
-void Force_Stress_LCAO<double>::calForceStressIntegralPart(const bool isGammaOnly,
+void Force_Stress_LCAO<double>::calForceStressIntegralPart(
+    const bool isGammaOnly,
     const bool isforce,
     const bool isstress,
     Local_Orbital_Charge& loc,
@@ -731,9 +738,12 @@ void Force_Stress_LCAO<double>::calForceStressIntegralPart(const bool isGammaOnl
 #else
     ModuleBase::matrix& svl_dphi,
 #endif
-    LCAO_Hamilt& uhm,
+    LCAO_Hamilt &uhm,
+	Parallel_Orbitals &pv,
+    LCAO_Matrix &lm,
     const K_Vectors& kv)
 {
+
     flk.ftable_gamma(isforce,
         isstress,
         psi,
@@ -752,11 +762,15 @@ void Force_Stress_LCAO<double>::calForceStressIntegralPart(const bool isGammaOnl
 #else
         svl_dphi,
 #endif
-        uhm);
+        uhm,
+        lm);
     return;
 }
+
+
 template<>
-void Force_Stress_LCAO<std::complex<double>>::calForceStressIntegralPart(const bool isGammaOnly,
+void Force_Stress_LCAO<std::complex<double>>::calForceStressIntegralPart(
+    const bool isGammaOnly,
     const bool isforce,
     const bool isstress,
     Local_Orbital_Charge& loc,
@@ -775,7 +789,9 @@ void Force_Stress_LCAO<std::complex<double>>::calForceStressIntegralPart(const b
 #else
     ModuleBase::matrix& svl_dphi,
 #endif
-    LCAO_Hamilt& uhm,
+	LCAO_Hamilt &uhm,
+	Parallel_Orbitals &pv,
+	LCAO_Matrix &lm,
     const K_Vectors& kv)
 {
         flk.ftable_k(isforce,
@@ -797,8 +813,10 @@ void Force_Stress_LCAO<std::complex<double>>::calForceStressIntegralPart(const b
 #else
                      svl_dphi,
 #endif
-                     uhm,
-            kv);
+					 uhm,
+					 pv,
+					 lm,
+                     kv);
     return;
 }
 
@@ -904,4 +922,4 @@ void Force_Stress_LCAO<T>::forceSymmetry(ModuleBase::matrix& fcs, ModuleSymmetry
 }
 
 template class Force_Stress_LCAO<double>;
-template class Force_Stress_LCAO<std::complex<double>>;
\ No newline at end of file
+template class Force_Stress_LCAO<std::complex<double>>;
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
index ccca8073ff6..8aaf0d4a387 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
@@ -30,8 +30,10 @@ class Force_Stress_LCAO
         const bool isstress,
         const bool istestf,
         const bool istests,
-        Local_Orbital_Charge& loc,
-        const elecstate::ElecState* pelec,
+		Local_Orbital_Charge& loc,
+		Parallel_Orbitals &pv,
+		LCAO_Matrix &lm,
+		const elecstate::ElecState* pelec,
         const psi::Psi<T>* psi,
         LCAO_Hamilt& uhm,
         ModuleBase::matrix& fcs,
@@ -85,8 +87,10 @@ class Force_Stress_LCAO
 #else
         ModuleBase::matrix& svl_dphi,
 #endif
-        LCAO_Hamilt& uhm,
-        const K_Vectors& kv);
+		LCAO_Hamilt &uhm,
+	    Parallel_Orbitals &pv,
+		LCAO_Matrix &lm,
+		const K_Vectors& kv);
 
     void calStressPwPart(ModuleBase::matrix& sigmadvl,
                          ModuleBase::matrix& sigmahar,
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
index 1db8883136d..2bc04ebd51b 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
@@ -38,7 +38,8 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce,
 #else
                                     ModuleBase::matrix& svl_dphi,
 #endif
-                                    LCAO_Hamilt& uhm)
+									LCAO_Hamilt &uhm,
+									LCAO_Matrix &lm)
 {
     ModuleBase::TITLE("Force_LCAO_gamma", "ftable");
     ModuleBase::timer::tick("Force_LCAO_gamma", "ftable_gamma");
@@ -53,15 +54,15 @@ 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);
+    this->allocate_gamma(*this->ParaV, lm);
 
     // calculate the 'energy density matrix' here.
-    this->cal_foverlap(isforce, isstress, psid, pelec, foverlap, soverlap);
+    this->cal_foverlap(isforce, isstress, psid, pelec, lm, foverlap, soverlap);
 
     // sum up the density matrix with different spin
     // DM->sum_DMR_spin();
     //
-    this->cal_ftvnl_dphi(DM, isforce, isstress, ftvnl_dphi, stvnl_dphi);
+    this->cal_ftvnl_dphi(DM, lm, isforce, isstress, ftvnl_dphi, stvnl_dphi);
     this->cal_fvnl_dbeta(DM, isforce, isstress, fvnl_dbeta, svnl_dbeta);
 
     this->cal_fvl_dphi(loc.DM, isforce, isstress, pelec->pot, fvl_dphi, svl_dphi);
@@ -119,13 +120,15 @@ void Force_LCAO_gamma::ftable_gamma(const bool isforce,
 
     // delete DSloc_x, DSloc_y, DSloc_z
     // delete DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z
-    this->finish_ftable_gamma();
+    this->finish_ftable_gamma(lm);
 
     ModuleBase::timer::tick("Force_LCAO_gamma", "ftable_gamma");
     return;
 }
 
-void Force_LCAO_gamma::allocate_gamma(const Parallel_Orbitals& pv)
+void Force_LCAO_gamma::allocate_gamma(
+    const Parallel_Orbitals &pv,
+    LCAO_Matrix &lm)
 {
     ModuleBase::TITLE("Force_LCAO_gamma", "allocate_gamma");
     ModuleBase::timer::tick("Force_LCAO_gamma", "allocate_gamma");
@@ -138,108 +141,134 @@ void Force_LCAO_gamma::allocate_gamma(const Parallel_Orbitals& pv)
     // liaochen add on 2010/7/12
     // save the results in dense matrix by now.
     // pv.nloc: number of H elements in this proc.
-    this->UHM->LM->DSloc_x = new double[pv.nloc];
-    this->UHM->LM->DSloc_y = new double[pv.nloc];
-    this->UHM->LM->DSloc_z = new double[pv.nloc];
-    ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_x, pv.nloc);
-    ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_y, pv.nloc);
-    ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_z, pv.nloc);
+    lm.DSloc_x = new double[pv.nloc];
+    lm.DSloc_y = new double[pv.nloc];
+    lm.DSloc_z = new double[pv.nloc];
+    ModuleBase::GlobalFunc::ZEROS(lm.DSloc_x, pv.nloc);
+    ModuleBase::GlobalFunc::ZEROS(lm.DSloc_y, pv.nloc);
+    ModuleBase::GlobalFunc::ZEROS(lm.DSloc_z, pv.nloc);
     ModuleBase::Memory::record("Force::dS_GO", sizeof(double) * pv.nloc * 3);
     // allocate stress part in gamma_only-line, added by zhengdy-stress
     if (GlobalV::CAL_STRESS)
     {
-        this->UHM->LM->DSloc_11 = new double[pv.nloc];
-        this->UHM->LM->DSloc_12 = new double[pv.nloc];
-        this->UHM->LM->DSloc_13 = new double[pv.nloc];
-        this->UHM->LM->DSloc_22 = new double[pv.nloc];
-        this->UHM->LM->DSloc_23 = new double[pv.nloc];
-        this->UHM->LM->DSloc_33 = new double[pv.nloc];
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_11, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_12, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_13, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_22, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_23, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_33, pv.nloc);
-        this->UHM->LM->DHloc_fixed_11 = new double[pv.nloc];
-        this->UHM->LM->DHloc_fixed_12 = new double[pv.nloc];
-        this->UHM->LM->DHloc_fixed_13 = new double[pv.nloc];
-        this->UHM->LM->DHloc_fixed_22 = new double[pv.nloc];
-        this->UHM->LM->DHloc_fixed_23 = new double[pv.nloc];
-        this->UHM->LM->DHloc_fixed_33 = new double[pv.nloc];
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_11, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_12, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_13, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_22, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_23, pv.nloc);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_33, pv.nloc);
+        lm.DSloc_11 = new double[pv.nloc];
+        lm.DSloc_12 = new double[pv.nloc];
+        lm.DSloc_13 = new double[pv.nloc];
+        lm.DSloc_22 = new double[pv.nloc];
+        lm.DSloc_23 = new double[pv.nloc];
+        lm.DSloc_33 = new double[pv.nloc];
+        ModuleBase::GlobalFunc::ZEROS(lm.DSloc_11, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DSloc_12, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DSloc_13, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DSloc_22, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DSloc_23, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DSloc_33, pv.nloc);
+        lm.DHloc_fixed_11 = new double[pv.nloc];
+        lm.DHloc_fixed_12 = new double[pv.nloc];
+        lm.DHloc_fixed_13 = new double[pv.nloc];
+        lm.DHloc_fixed_22 = new double[pv.nloc];
+        lm.DHloc_fixed_23 = new double[pv.nloc];
+        lm.DHloc_fixed_33 = new double[pv.nloc];
+        ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixed_11, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixed_12, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixed_13, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixed_22, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixed_23, pv.nloc);
+        ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixed_33, pv.nloc);
         ModuleBase::Memory::record("Stress::dSH_GO", sizeof(double) * pv.nloc * 12);
     }
     // calculate dS in LCAO basis
     // ModuleBase::timer::tick("Force_LCAO_gamma","build_S_new");
-    this->UHM->genH.build_ST_new('S', cal_deri, GlobalC::ucell, this->UHM->LM->Sloc.data());
+    this->UHM->genH.build_ST_new('S', cal_deri, GlobalC::ucell, lm.Sloc.data());
     // ModuleBase::timer::tick("Force_LCAO_gamma","build_S_new");
 
     // calculate dT in LCAP
     // allocation dt
     // liaochen add on 2010/7/12
-    this->UHM->LM->DHloc_fixed_x = new double[pv.nloc];
-    this->UHM->LM->DHloc_fixed_y = new double[pv.nloc];
-    this->UHM->LM->DHloc_fixed_z = new double[pv.nloc];
+    lm.DHloc_fixed_x = new double[pv.nloc];
+    lm.DHloc_fixed_y = new double[pv.nloc];
+    lm.DHloc_fixed_z = new double[pv.nloc];
     ModuleBase::Memory::record("Force::dTVNL", sizeof(double) * pv.nloc * 3);
-    ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_x, pv.nloc);
-    ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_y, pv.nloc);
-    ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_z, pv.nloc);
+    ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixed_x, pv.nloc);
+    ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixed_y, pv.nloc);
+    ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixed_z, pv.nloc);
 
     // calculate dT
     // calculate T + VNL(P1) in LCAO basis
     // ModuleBase::timer::tick("Force_LCAO_gamma","build_T_new");
-    this->UHM->genH.build_ST_new('T', cal_deri, GlobalC::ucell, this->UHM->LM->Hloc_fixed.data());
+    this->UHM->genH.build_ST_new('T', cal_deri, GlobalC::ucell, lm.Hloc_fixed.data());
     // ModuleBase::timer::tick("Force_LCAO_gamma","build_T_new");
-    // test_gamma(this->UHM->LM->DHloc_fixed_x, "dHloc_fixed_x T part");
+    // test_gamma(lm.DHloc_fixed_x, "dHloc_fixed_x T part");
 
     // ModuleBase::timer::tick("Force_LCAO_gamma","build_Nonlocal_mu");
-    this->UHM->genH.build_Nonlocal_mu_new(this->UHM->genH.LM->Hloc_fixed.data(), cal_deri);
+    this->UHM->genH.build_Nonlocal_mu_new(lm.Hloc_fixed.data(), cal_deri);
     // ModuleBase::timer::tick("Force_LCAO_gamma","build_Nonlocal_mu");
-    // test_gamma(this->UHM->LM->DHloc_fixed_x, "dHloc_fixed_x Vnl part");
+    // test_gamma(lm.DHloc_fixed_x, "dHloc_fixed_x Vnl part");
 
     // calculate asynchronous S matrix to output for Hefei-NAMD
     if (INPUT.cal_syns)
     {
         cal_deri = false;
-        this->UHM->genH.LM->zeros_HSgamma('S');
-        this->UHM->genH
-            .build_ST_new('S', cal_deri, GlobalC::ucell, this->UHM->genH.LM->Sloc.data(), INPUT.cal_syns, INPUT.dmax);
+
+        lm.zeros_HSgamma('S');
+
+        this->UHM->genH.build_ST_new('S', cal_deri, GlobalC::ucell, 
+
+        lm.Sloc.data(), INPUT.cal_syns, INPUT.dmax);
+
         bool bit = false; // LiuXh, 2017-03-21
-        ModuleIO::save_mat(0, this->UHM->genH.LM->Hloc.data(), GlobalV::NLOCAL, bit, GlobalV::out_ndigits, 0, GlobalV::out_app_flag, "H", "data-" + std::to_string(0), *this->ParaV, GlobalV::DRANK);
-        ModuleIO::save_mat(0, this->UHM->genH.LM->Sloc.data(), GlobalV::NLOCAL, bit, GlobalV::out_ndigits, 0, GlobalV::out_app_flag, "S", "data-" + std::to_string(0), *this->ParaV, GlobalV::DRANK);
-    }
+
+		ModuleIO::save_mat(0, 
+				lm.Hloc.data(), 
+				GlobalV::NLOCAL, 
+				bit, 
+				GlobalV::out_ndigits, 
+				0, 
+				GlobalV::out_app_flag, 
+				"H", 
+				"data-" + std::to_string(0), 
+				*this->ParaV, 
+				GlobalV::DRANK);
+
+		ModuleIO::save_mat(0, 
+				lm.Sloc.data(), 
+				GlobalV::NLOCAL, 
+				bit, 
+				GlobalV::out_ndigits, 
+				0, 
+				GlobalV::out_app_flag, 
+				"S", 
+				"data-" + std::to_string(0), 
+				*this->ParaV, 
+				GlobalV::DRANK);
+	}
 
     ModuleBase::timer::tick("Force_LCAO_gamma", "allocate_gamma");
     return;
 }
 
-void Force_LCAO_gamma::finish_ftable_gamma(void)
+void Force_LCAO_gamma::finish_ftable_gamma(LCAO_Matrix &lm)
 {
-    delete[] this->UHM->LM->DSloc_x;
-    delete[] this->UHM->LM->DSloc_y;
-    delete[] this->UHM->LM->DSloc_z;
-    delete[] this->UHM->LM->DHloc_fixed_x;
-    delete[] this->UHM->LM->DHloc_fixed_y;
-    delete[] this->UHM->LM->DHloc_fixed_z;
+    delete[] lm.DSloc_x;
+    delete[] lm.DSloc_y;
+    delete[] lm.DSloc_z;
+    delete[] lm.DHloc_fixed_x;
+    delete[] lm.DHloc_fixed_y;
+    delete[] lm.DHloc_fixed_z;
     if (GlobalV::CAL_STRESS) // added by zhengdy-stress
     {
-        delete[] this->UHM->LM->DSloc_11;
-        delete[] this->UHM->LM->DSloc_12;
-        delete[] this->UHM->LM->DSloc_13;
-        delete[] this->UHM->LM->DHloc_fixed_11;
-        delete[] this->UHM->LM->DHloc_fixed_12;
-        delete[] this->UHM->LM->DHloc_fixed_13;
-        delete[] this->UHM->LM->DSloc_22;
-        delete[] this->UHM->LM->DSloc_23;
-        delete[] this->UHM->LM->DSloc_33;
-        delete[] this->UHM->LM->DHloc_fixed_22;
-        delete[] this->UHM->LM->DHloc_fixed_23;
-        delete[] this->UHM->LM->DHloc_fixed_33;
+        delete[] lm.DSloc_11;
+        delete[] lm.DSloc_12;
+        delete[] lm.DSloc_13;
+        delete[] lm.DHloc_fixed_11;
+        delete[] lm.DHloc_fixed_12;
+        delete[] lm.DHloc_fixed_13;
+        delete[] lm.DSloc_22;
+        delete[] lm.DSloc_23;
+        delete[] lm.DSloc_33;
+        delete[] lm.DHloc_fixed_22;
+        delete[] lm.DHloc_fixed_23;
+        delete[] lm.DHloc_fixed_33;
     }
     return;
 }
@@ -282,4 +311,4 @@ void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix&
         }
     }
 }
-} // namespace StressTools
\ No newline at end of file
+} // namespace StressTools
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h
index e9efd88b425..a3358382ee7 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.h
@@ -20,8 +20,11 @@ class Force_LCAO_gamma
     ~Force_LCAO_gamma();
 
   private:
+
     LCAO_Hamilt* UHM;
+
     const Parallel_Orbitals* ParaV;
+
     elecstate::Potential* pot;
 
     // orthonormal force + contribution from T and VNL
@@ -43,12 +46,17 @@ class Force_LCAO_gamma
 #else
                       ModuleBase::matrix& svl_dphi,
 #endif
-                      LCAO_Hamilt& uhm);
+                      LCAO_Hamilt& uhm,
+                      LCAO_Matrix& lm);
+
+
 
     // get the ds, dt, dvnl.
-    void allocate_gamma(const Parallel_Orbitals& pv);
+	void allocate_gamma(
+			const Parallel_Orbitals& pv,
+			LCAO_Matrix &lm);
 
-    void finish_ftable_gamma(void);
+    void finish_ftable_gamma(LCAO_Matrix &lm);
 
     void average_force(double* fm);
 
@@ -61,15 +69,17 @@ class Force_LCAO_gamma
 
     void cal_foverlap(const bool isforce,
                       const bool isstress,
-                      const psi::Psi<double>* psid,
-                      const elecstate::ElecState* pelec,
-                      ModuleBase::matrix& foverlap,
-                      ModuleBase::matrix& soverlap);
+                      const psi::Psi<double> *psid,
+                      const elecstate::ElecState *pelec,
+                      LCAO_Matrix &lm,
+                      ModuleBase::matrix &foverlap,
+                      ModuleBase::matrix &soverlap);
 
     //-------------------------------------------------------------
     // forces related to kinetic and non-local pseudopotentials
     //--------------------------------------------------------------
     void cal_ftvnl_dphi(const elecstate::DensityMatrix<double, double>* DM,
+                        LCAO_Matrix &lm,
                         const bool isforce,
                         const bool isstress,
                         ModuleBase::matrix& ftvnl_dphi,
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_edm.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_edm.cpp
index a69821c50ce..e719bc8537e 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_edm.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_edm.cpp
@@ -11,10 +11,11 @@
 void Force_LCAO_gamma::cal_foverlap(
 	const bool isforce, 
     const bool isstress,
-    const psi::Psi<double>* psid,
-    const elecstate::ElecState* pelec,
-    ModuleBase::matrix& foverlap,
-	ModuleBase::matrix& soverlap)
+    const psi::Psi<double> *psid,
+    const elecstate::ElecState *pelec,
+    LCAO_Matrix &lm,
+    ModuleBase::matrix &foverlap,
+	ModuleBase::matrix &soverlap)
 {
     ModuleBase::TITLE("Force_LCAO_gamma","cal_foverlap");
     ModuleBase::timer::tick("Force_LCAO_gamma","cal_foverlap");
@@ -61,19 +62,19 @@ void Force_LCAO_gamma::cal_foverlap(
 
                 if(isforce)
                 {
-                    foverlap(iat,0) += sum * this->UHM->LM->DSloc_x[index];
-                    foverlap(iat,1) += sum * this->UHM->LM->DSloc_y[index];
-                    foverlap(iat,2) += sum * this->UHM->LM->DSloc_z[index];
+                    foverlap(iat,0) += sum * lm.DSloc_x[index];
+                    foverlap(iat,1) += sum * lm.DSloc_y[index];
+                    foverlap(iat,2) += sum * lm.DSloc_z[index];
                 }
 
                 if(isstress)
                 {
-                    soverlap(0,0) += sum/2.0 * this->UHM->LM->DSloc_11[index];
-                    soverlap(0,1) += sum/2.0 * this->UHM->LM->DSloc_12[index];
-                    soverlap(0,2) += sum/2.0 * this->UHM->LM->DSloc_13[index];
-                    soverlap(1,1) += sum/2.0 * this->UHM->LM->DSloc_22[index];
-                    soverlap(1,2) += sum/2.0 * this->UHM->LM->DSloc_23[index];
-                    soverlap(2,2) += sum/2.0 * this->UHM->LM->DSloc_33[index];
+                    soverlap(0,0) += sum/2.0 * lm.DSloc_11[index];
+                    soverlap(0,1) += sum/2.0 * lm.DSloc_12[index];
+                    soverlap(0,2) += sum/2.0 * lm.DSloc_13[index];
+                    soverlap(1,1) += sum/2.0 * lm.DSloc_22[index];
+                    soverlap(1,2) += sum/2.0 * lm.DSloc_23[index];
+                    soverlap(2,2) += sum/2.0 * lm.DSloc_33[index];
                 }
             }
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_tvnl.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_tvnl.cpp
index 3759f5f9a7e..dd89912d1dd 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_tvnl.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma_tvnl.cpp
@@ -237,6 +237,7 @@ void Force_LCAO_gamma::cal_fvnl_dbeta(
 
 void Force_LCAO_gamma::cal_ftvnl_dphi(
     const elecstate::DensityMatrix<double, double>* DM,
+    LCAO_Matrix &lm,
     const bool isforce, 
 	const bool isstress, 
 	ModuleBase::matrix& ftvnl_dphi, 
@@ -268,25 +269,26 @@ void Force_LCAO_gamma::cal_ftvnl_dphi(
 
                 if(isforce)
 				{
-					ftvnl_dphi(iat,0) += sum * this->UHM->LM->DHloc_fixed_x[index];
-					ftvnl_dphi(iat,1) += sum * this->UHM->LM->DHloc_fixed_y[index];
-					ftvnl_dphi(iat,2) += sum * this->UHM->LM->DHloc_fixed_z[index];
+					ftvnl_dphi(iat,0) += sum * lm.DHloc_fixed_x[index];
+					ftvnl_dphi(iat,1) += sum * lm.DHloc_fixed_y[index];
+					ftvnl_dphi(iat,2) += sum * lm.DHloc_fixed_z[index];
 				}
                 if(isstress)
                 {
-                    stvnl_dphi(0,0) += sum/2.0 * this->UHM->LM->DHloc_fixed_11[index];
-                    stvnl_dphi(0,1) += sum/2.0 * this->UHM->LM->DHloc_fixed_12[index];
-                    stvnl_dphi(0,2) += sum/2.0 * this->UHM->LM->DHloc_fixed_13[index];
-                    stvnl_dphi(1,1) += sum/2.0 * this->UHM->LM->DHloc_fixed_22[index];
-                    stvnl_dphi(1,2) += sum/2.0 * this->UHM->LM->DHloc_fixed_23[index];
-                    stvnl_dphi(2,2) += sum/2.0 * this->UHM->LM->DHloc_fixed_33[index];   
+                    stvnl_dphi(0,0) += sum/2.0 * lm.DHloc_fixed_11[index];
+                    stvnl_dphi(0,1) += sum/2.0 * lm.DHloc_fixed_12[index];
+                    stvnl_dphi(0,2) += sum/2.0 * lm.DHloc_fixed_13[index];
+                    stvnl_dphi(1,1) += sum/2.0 * lm.DHloc_fixed_22[index];
+                    stvnl_dphi(1,2) += sum/2.0 * lm.DHloc_fixed_23[index];
+                    stvnl_dphi(2,2) += sum/2.0 * lm.DHloc_fixed_33[index];   
                 }
             }
         }
     }
-    if(isstress){
-        StressTools::stress_fill(GlobalC::ucell.lat0, GlobalC::ucell.omega, stvnl_dphi);
-    }
+	if(isstress)
+	{
+		StressTools::stress_fill(GlobalC::ucell.lat0, GlobalC::ucell.omega, stvnl_dphi);
+	}
     ModuleBase::timer::tick("Force_LCAO_gamma","cal_ftvnl_dphi");
     return;
-}
\ No newline at end of file
+}
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
index 918d4c53461..947d5bc5be8 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
@@ -50,7 +50,9 @@ void Force_LCAO_k::ftable_k(const bool isforce,
 #else
                             ModuleBase::matrix& svl_dphi,
 #endif
-                            LCAO_Hamilt& uhm,
+                            LCAO_Hamilt &uhm,
+                            Parallel_Orbitals &pv,
+                            LCAO_Matrix &lm,
                             const K_Vectors& kv)
 {
     ModuleBase::TITLE("Force_LCAO_k", "ftable_k");
@@ -61,25 +63,22 @@ void Force_LCAO_k::ftable_k(const bool isforce,
     elecstate::DensityMatrix<complex<double>,double>* DM
     = dynamic_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(pelec)->get_DM();
 
-    this->ParaV = DM->get_paraV_pointer();
-    const Parallel_Orbitals* pv = this->ParaV;
-    //const Parallel_Orbitals* pv = loc.ParaV;
-    this->allocate_k(*this->ParaV, kv.nks, kv.kvec_d);
+    this->allocate_k(pv, lm, kv.nks, kv.kvec_d);
 
     // calculate the energy density matrix
     // and the force related to overlap matrix and energy density matrix.
-    this->cal_foverlap_k(isforce, isstress, ra, psi, loc, DM, foverlap, soverlap, pelec, kv.nks, kv);
+    this->cal_foverlap_k(isforce, isstress, ra, psi, loc, pv, lm, DM, foverlap, soverlap, pelec, kv.nks, kv);
 
     //DM->sum_DMR_spin();
 
-    this->cal_ftvnl_dphi_k(DM, isforce, isstress, ra, ftvnl_dphi, stvnl_dphi);
+    this->cal_ftvnl_dphi_k(DM, pv, lm, isforce, isstress, ra, ftvnl_dphi, stvnl_dphi);
 
     // ---------------------------------------
     // doing on the real space grid.
     // ---------------------------------------
-    this->cal_fvl_dphi_k(isforce, isstress, pelec->pot, fvl_dphi, svl_dphi, loc.DM_R);
+    this->cal_fvl_dphi_k(isforce, isstress, lm, pelec->pot, fvl_dphi, svl_dphi, loc.DM_R);
 
-    this->cal_fvnl_dbeta_k(DM, isforce, isstress, fvnl_dbeta, svnl_dbeta);
+    this->cal_fvnl_dbeta_k(DM, isforce, isstress, pv, fvnl_dbeta, svnl_dbeta);
 
 #ifdef __DEEPKS
     if (GlobalV::deepks_scf)
@@ -149,49 +148,58 @@ void Force_LCAO_k::ftable_k(const bool isforce,
 }
 
 void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv,
+                              LCAO_Matrix &lm,
                               const int& nks,
                               const std::vector<ModuleBase::Vector3<double>>& kvec_d)
 {
     ModuleBase::TITLE("Force_LCAO_k", "allocate_k");
     ModuleBase::timer::tick("Force_LCAO_k", "allocate_k");
 
-    this->ParaV = &pv;
     const int nnr = pv.nnr;
+
+    assert(nnr>=0);
+
     //--------------------------------
     // (1) allocate for dSx dSy & dSz
     //--------------------------------
-    this->UHM->LM->DSloc_Rx = new double[nnr];
-    this->UHM->LM->DSloc_Ry = new double[nnr];
-    this->UHM->LM->DSloc_Rz = new double[nnr];
-    const auto init_DSloc_Rxyz = [this, 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(this->UHM->LM->DSloc_Rx + beg, len);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_Ry + beg, len);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DSloc_Rz + beg, len);
-    };
+    lm.DSloc_Rx = new double[nnr];
+    lm.DSloc_Ry = new double[nnr];
+    lm.DSloc_Rz = new double[nnr];
+
+    // mohan add lm on 2024-03-31
+	const auto init_DSloc_Rxyz = [this, nnr, &lm](int num_threads, int thread_id) 
+	{
+		int beg=0;
+		int len=0;
+		ModuleBase::BLOCK_TASK_DIST_1D(num_threads, thread_id, nnr, 1024, 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);
     ModuleBase::Memory::record("Force::dS_K", sizeof(double) * nnr * 3);
 
     if (GlobalV::CAL_STRESS)
     {
-        this->UHM->LM->DH_r = new double[3 * nnr];
-        this->UHM->LM->stvnl11 = new double[nnr];
-        this->UHM->LM->stvnl12 = new double[nnr];
-        this->UHM->LM->stvnl13 = new double[nnr];
-        this->UHM->LM->stvnl22 = new double[nnr];
-        this->UHM->LM->stvnl23 = new double[nnr];
-        this->UHM->LM->stvnl33 = new double[nnr];
-        const auto init_DH_r_stvnl = [this, nnr](int num_threads, int thread_id) {
+        lm.DH_r = new double[3 * nnr];
+        lm.stvnl11 = new double[nnr];
+        lm.stvnl12 = new double[nnr];
+        lm.stvnl13 = new double[nnr];
+        lm.stvnl22 = new double[nnr];
+        lm.stvnl23 = new double[nnr];
+        lm.stvnl33 = new double[nnr];
+        // mohan add lm on 2024-03-31
+        const auto init_DH_r_stvnl = [this, nnr, &lm](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(this->UHM->LM->DH_r + 3 * beg, 3 * len);
-            ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->stvnl11 + beg, len);
-            ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->stvnl12 + beg, len);
-            ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->stvnl13 + beg, len);
-            ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->stvnl22 + beg, len);
-            ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->stvnl23 + beg, len);
-            ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->stvnl33 + beg, len);
+            ModuleBase::GlobalFunc::ZEROS(lm.DH_r + 3 * beg, 3 * len);
+            ModuleBase::GlobalFunc::ZEROS(lm.stvnl11 + beg, len);
+            ModuleBase::GlobalFunc::ZEROS(lm.stvnl12 + beg, len);
+            ModuleBase::GlobalFunc::ZEROS(lm.stvnl13 + beg, len);
+            ModuleBase::GlobalFunc::ZEROS(lm.stvnl22 + beg, len);
+            ModuleBase::GlobalFunc::ZEROS(lm.stvnl23 + beg, len);
+            ModuleBase::GlobalFunc::ZEROS(lm.stvnl33 + beg, len);
         };
         ModuleBase::OMP_PARALLEL(init_DH_r_stvnl);
 
@@ -208,15 +216,17 @@ void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv,
     //-----------------------------------------
     // (2) allocate for <phi | T + Vnl | dphi>
     //-----------------------------------------
-    this->UHM->LM->DHloc_fixedR_x = new double[nnr];
-    this->UHM->LM->DHloc_fixedR_y = new double[nnr];
-    this->UHM->LM->DHloc_fixedR_z = new double[nnr];
-    const auto init_DHloc_fixedR_xyz = [this, nnr](int num_threads, int thread_id) {
+    lm.DHloc_fixedR_x = new double[nnr];
+    lm.DHloc_fixedR_y = new double[nnr];
+    lm.DHloc_fixedR_z = new double[nnr];
+
+    // mohan add lm on 2024-03-31
+    const auto init_DHloc_fixedR_xyz = [this, nnr, &lm](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(this->UHM->LM->DHloc_fixedR_x + beg, len);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixedR_y + beg, len);
-        ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixedR_z + beg, len);
+        ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixedR_x + beg, len);
+        ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixedR_y + beg, len);
+        ModuleBase::GlobalFunc::ZEROS(lm.DHloc_fixedR_z + beg, len);
     };
     ModuleBase::OMP_PARALLEL(init_DHloc_fixedR_xyz);
     ModuleBase::Memory::record("Force::dTVNL", sizeof(double) * nnr * 3);
@@ -224,27 +234,51 @@ void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv,
     // calculate dT=<phi|kin|dphi> in LCAO
     // calculate T + VNL(P1) in LCAO basis
     this->UHM->genH.build_ST_new('T', cal_deri, GlobalC::ucell, this->UHM->genH.LM->Hloc_fixedR.data());
-    // test(this->UHM->LM->DHloc_fixedR_x,"this->UHM->LM->DHloc_fixedR_x T part");
 
     // calculate dVnl=<phi|dVnl|dphi> in LCAO
     this->UHM->genH.build_Nonlocal_mu_new(this->UHM->genH.LM->Hloc_fixed.data(), cal_deri);
-    // test(this->UHM->LM->DHloc_fixedR_x,"this->UHM->LM->DHloc_fixedR_x Vnl part");
 
     // calculate asynchronous S matrix to output for Hefei-NAMD
     if (INPUT.cal_syns)
     {
         cal_deri = false;
         // this->UHM->genH.build_ST_new('S', cal_deri, GlobalC::ucell, this->UHM->genH.LM->SlocR.data(),
-        // INPUT.cal_syns);
-        this->UHM->genH
-            .build_ST_new('S', cal_deri, GlobalC::ucell, this->UHM->genH.LM->SlocR.data(), INPUT.cal_syns, INPUT.dmax);
+		// INPUT.cal_syns);
+		this->UHM->genH.build_ST_new('S', 
+				cal_deri, 
+				GlobalC::ucell, 
+				lm.SlocR.data(), 
+				INPUT.cal_syns, 
+				INPUT.dmax);
+
         for (int ik = 0; ik < nks; ik++)
         {
-            this->UHM->genH.LM->zeros_HSk('S');
-            this->UHM->genH.LM->folding_fixedH(ik, kvec_d, 1);
+            lm.zeros_HSk('S');
+            lm.folding_fixedH(ik, kvec_d, 1);
             bool bit = false; // LiuXh, 2017-03-21
-            ModuleIO::save_mat(0, this->UHM->genH.LM->Hloc2.data(), GlobalV::NLOCAL, bit, GlobalV::out_ndigits, 0, GlobalV::out_app_flag, "H", "data-" + std::to_string(ik), *this->ParaV, GlobalV::DRANK);
-            ModuleIO::save_mat(0, this->UHM->genH.LM->Sloc2.data(), GlobalV::NLOCAL, bit, GlobalV::out_ndigits, 0, GlobalV::out_app_flag, "S", "data-" + std::to_string(ik), *this->ParaV, GlobalV::DRANK);
+			ModuleIO::save_mat(0, 
+					lm.Hloc2.data(), 
+					GlobalV::NLOCAL, 
+					bit, 
+					GlobalV::out_ndigits, 
+					0, 
+					GlobalV::out_app_flag, 
+					"H", 
+					"data-" + std::to_string(ik), 
+					pv, 
+					GlobalV::DRANK);
+
+            ModuleIO::save_mat(0, 
+					lm.Sloc2.data(), 
+					GlobalV::NLOCAL, 
+					bit, 
+					GlobalV::out_ndigits, 
+					0, 
+					GlobalV::out_app_flag, 
+					"S", 
+					"data-" + std::to_string(ik), 
+					pv, 
+					GlobalV::DRANK);
         }
     }
 
@@ -252,23 +286,25 @@ void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv,
     return;
 }
 
-void Force_LCAO_k::finish_k(void)
+
+void Force_LCAO_k::finish_k(LCAO_Matrix &lm)
 {
-    delete[] this->UHM->LM->DSloc_Rx;
-    delete[] this->UHM->LM->DSloc_Ry;
-    delete[] this->UHM->LM->DSloc_Rz;
-    delete[] this->UHM->LM->DHloc_fixedR_x;
-    delete[] this->UHM->LM->DHloc_fixedR_y;
-    delete[] this->UHM->LM->DHloc_fixedR_z;
+    delete[] lm.DSloc_Rx;
+    delete[] lm.DSloc_Ry;
+    delete[] lm.DSloc_Rz;
+    delete[] lm.DHloc_fixedR_x;
+    delete[] lm.DHloc_fixedR_y;
+    delete[] lm.DHloc_fixedR_z;
+
     if (GlobalV::CAL_STRESS)
     {
-        delete[] this->UHM->LM->DH_r;
-        delete[] this->UHM->LM->stvnl11;
-        delete[] this->UHM->LM->stvnl12;
-        delete[] this->UHM->LM->stvnl13;
-        delete[] this->UHM->LM->stvnl22;
-        delete[] this->UHM->LM->stvnl23;
-        delete[] this->UHM->LM->stvnl33;
+        delete[] lm.DH_r;
+        delete[] lm.stvnl11;
+        delete[] lm.stvnl12;
+        delete[] lm.stvnl13;
+        delete[] lm.stvnl22;
+        delete[] lm.stvnl23;
+        delete[] lm.stvnl33;
     }
     return;
 }
@@ -276,23 +312,23 @@ void Force_LCAO_k::finish_k(void)
 #include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h"
 void Force_LCAO_k::cal_foverlap_k(const bool isforce,
                                   const bool isstress,
-                                  Record_adj& ra,
-                                  const psi::Psi<std::complex<double>>* psi,
-                                  Local_Orbital_Charge& loc,
-                                  const elecstate::DensityMatrix<std::complex<double>, double>* DM,
-                                  ModuleBase::matrix& foverlap,
-                                  ModuleBase::matrix& soverlap,
-                                  const elecstate::ElecState* pelec,
-                                  const int& nks,
-                                  const K_Vectors& kv)
+                                  Record_adj &ra,
+                                  const psi::Psi<std::complex<double>> *psi,
+                                  Local_Orbital_Charge &loc,
+                                  Parallel_Orbitals &pv,
+                                  LCAO_Matrix &lm,
+                                  const elecstate::DensityMatrix<std::complex<double>, double> *DM,
+                                  ModuleBase::matrix &foverlap,
+                                  ModuleBase::matrix &soverlap,
+                                  const elecstate::ElecState *pelec,
+                                  const int &nks,
+                                  const K_Vectors &kv)
 {
     ModuleBase::TITLE("Force_LCAO_k", "cal_foverlap_k");
     ModuleBase::timer::tick("Force_LCAO_k", "cal_foverlap_k");
 
-    const Parallel_Orbitals* pv = this->ParaV;
-
     // construct a DensityMatrix object
-    elecstate::DensityMatrix<std::complex<double>, double> EDM(&kv,pv,GlobalV::NSPIN);
+    elecstate::DensityMatrix<std::complex<double>, double> EDM(&kv,&pv,GlobalV::NSPIN);
     
     //--------------------------------------------
     // calculate the energy density matrix here.
@@ -328,7 +364,7 @@ void Force_LCAO_k::cal_foverlap_k(const bool isforce,
     }
     else
     {
-        //elecstate::cal_dm(loc.ParaV, wgEkb, psi[0], edm_k);
+        //elecstate::cal_dm(pv, wgEkb, psi[0], edm_k);
         // cal_dm_psi
         elecstate::cal_dm_psi(EDM.get_paraV_pointer(), wgEkb, psi[0], EDM);
     }
@@ -348,16 +384,16 @@ void Force_LCAO_k::cal_foverlap_k(const bool isforce,
     int total_irr = 0;
 #ifdef _OPENMP
 #pragma omp parallel
-    {
-        int num_threads = omp_get_num_threads();
-        ModuleBase::matrix local_soverlap(3, 3);
-        int local_total_irr = 0;
+	{
+		int num_threads = omp_get_num_threads();
+		ModuleBase::matrix local_soverlap(3, 3);
+		int local_total_irr = 0;
 #else
-    ModuleBase::matrix& local_soverlap = soverlap;
-    int& local_total_irr = total_irr;
+		ModuleBase::matrix& local_soverlap = soverlap;
+		int& local_total_irr = total_irr;
 #endif
 
-        ModuleBase::Vector3<double> tau1, dtau, tau2;
+		ModuleBase::Vector3<double> tau1, dtau, tau2;
 
 #ifdef _OPENMP
 #pragma omp for schedule(dynamic)
@@ -370,8 +406,10 @@ void Force_LCAO_k::cal_foverlap_k(const bool isforce,
             // get iat1
             int iat1 = GlobalC::ucell.itia2iat(T1, I1);
             double* foverlap_iat;
-            if (isforce)
-                foverlap_iat = &foverlap(iat, 0);
+			if (isforce)
+			{
+				foverlap_iat = &foverlap(iat, 0);
+			}
 
 #ifdef _OPENMP
             // using local stack to avoid false sharing in multi-threaded case
@@ -381,7 +419,7 @@ void Force_LCAO_k::cal_foverlap_k(const bool isforce,
                 foverlap_iat = foverlap_temp;
             }
 #endif
-            int irr = pv->nlocstart[iat];
+            int irr = pv.nlocstart[iat];
             const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0);
             for (int cb = 0; cb < ra.na_each[iat]; ++cb)
             {
@@ -399,13 +437,16 @@ void Force_LCAO_k::cal_foverlap_k(const bool isforce,
                 double Rz = ra.info[iat][cb][2];
                 // get BaseMatrix
                 hamilt::BaseMatrix<double>* tmp_matrix = EDM.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz);
-                if(tmp_matrix == nullptr) continue;
-                int row_ap = pv->atom_begin_row[iat1];
-                int col_ap = pv->atom_begin_col[iat2];
+				if(tmp_matrix == nullptr) 
+				{
+					continue;
+				}
+                int row_ap = pv.atom_begin_row[iat1];
+                int col_ap = pv.atom_begin_col[iat2];
                 // get DMR
-                for (int mu = 0; mu < pv->get_row_size(iat1); ++mu)
+                for (int mu = 0; mu < pv.get_row_size(iat1); ++mu)
                 {
-                    for (int nu = 0; nu < pv->get_col_size(iat2); ++nu)
+                    for (int nu = 0; nu < pv.get_col_size(iat2); ++nu)
                     {
                         // here do not sum over spin due to EDM.sum_DMR_spin();
                         double edm2d1 = tmp_matrix->get_value(mu,nu);
@@ -413,24 +454,28 @@ void Force_LCAO_k::cal_foverlap_k(const bool isforce,
 
                         if (isforce)
                         {
-                            foverlap_iat[0] -= edm2d2 * this->UHM->LM->DSloc_Rx[irr];
-                            foverlap_iat[1] -= edm2d2 * this->UHM->LM->DSloc_Ry[irr];
-                            foverlap_iat[2] -= edm2d2 * this->UHM->LM->DSloc_Rz[irr];
+                            foverlap_iat[0] -= edm2d2 * lm.DSloc_Rx[irr];
+                            foverlap_iat[1] -= edm2d2 * lm.DSloc_Ry[irr];
+                            foverlap_iat[2] -= edm2d2 * lm.DSloc_Rz[irr];
                         }
                         if (isstress)
                         {
                             for (int ipol = 0; ipol < 3; ipol++)
                             {
-                                local_soverlap(0, ipol) += edm2d1 * this->UHM->LM->DSloc_Rx[irr]
-                                                            * this->UHM->LM->DH_r[irr * 3 + ipol];
+                                local_soverlap(0, ipol) += edm2d1 * lm.DSloc_Rx[irr]
+                                                            * lm.DH_r[irr * 3 + ipol];
                                 if (ipol < 1)
-                                    continue;
-                                local_soverlap(1, ipol) += edm2d1 * this->UHM->LM->DSloc_Ry[irr]
-                                                            * this->UHM->LM->DH_r[irr * 3 + ipol];
-                                if (ipol < 2)
-                                    continue;
-                                local_soverlap(2, ipol) += edm2d1 * this->UHM->LM->DSloc_Rz[irr]
-                                                            * this->UHM->LM->DH_r[irr * 3 + ipol];
+								{
+									continue;
+								}
+                                local_soverlap(1, ipol) += edm2d1 * lm.DSloc_Ry[irr]
+                                                            * lm.DH_r[irr * 3 + ipol];
+								if (ipol < 2)
+								{
+									continue;
+								}
+                                local_soverlap(2, ipol) += edm2d1 * lm.DSloc_Rz[irr]
+                                                            * lm.DH_r[irr * 3 + ipol];
                             }
                         }
                         //}
@@ -457,12 +502,16 @@ void Force_LCAO_k::cal_foverlap_k(const bool isforce,
                 for (int ipol = 0; ipol < 3; ipol++)
                 {
                     soverlap(0, ipol) += local_soverlap(0, ipol);
-                    if (ipol < 1)
-                        continue;
-                    soverlap(1, ipol) += local_soverlap(1, ipol);
-                    if (ipol < 2)
-                        continue;
-                    soverlap(2, ipol) += local_soverlap(2, ipol);
+					if (ipol < 1)
+					{
+						continue;
+					}
+					soverlap(1, ipol) += local_soverlap(1, ipol);
+					if (ipol < 2)
+					{
+						continue;
+					}
+					soverlap(2, ipol) += local_soverlap(2, ipol);
                 }
             }
         }
@@ -474,10 +523,10 @@ void Force_LCAO_k::cal_foverlap_k(const bool isforce,
         StressTools::stress_fill(GlobalC::ucell.lat0, GlobalC::ucell.omega, soverlap);
     }
 
-    if (total_irr != pv->nnr)
+    if (total_irr != pv.nnr)
     {
         ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wrong irr", total_irr);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wrong LNNR.nnr", pv->nnr);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wrong LNNR.nnr", pv.nnr);
         ModuleBase::WARNING_QUIT("Force_LCAO_k::cal_foverlap_k", "irr!=LNNR.nnr");
     }
 
@@ -486,6 +535,8 @@ void Force_LCAO_k::cal_foverlap_k(const bool isforce,
 }
 
 void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<double>, double>* DM,
+                                    const Parallel_Orbitals &pv,
+                                    LCAO_Matrix &lm,
                                     const bool isforce,
                                     const bool isstress,
                                     Record_adj& ra,
@@ -496,7 +547,6 @@ void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<
     ModuleBase::timer::tick("Force_LCAO_k", "cal_ftvnl_dphi_k");
 
     int total_irr = 0;
-    const Parallel_Orbitals* pv = this->ParaV;
     // get the adjacent atom's information.
 
     //	GlobalV::ofs_running << " calculate the ftvnl_dphi_k force" << std::endl;
@@ -508,10 +558,10 @@ void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<
         int local_total_irr = 0;
 #pragma omp for schedule(dynamic)
 #else
-    ModuleBase::matrix& local_stvnl_dphi = stvnl_dphi;
-    int& local_total_irr = total_irr;
+		ModuleBase::matrix& local_stvnl_dphi = stvnl_dphi;
+		int& local_total_irr = total_irr;
 #endif
-        for (int iat = 0; iat < GlobalC::ucell.nat; iat++)
+		for (int iat = 0; iat < GlobalC::ucell.nat; iat++)
         {
             const int T1 = GlobalC::ucell.iat2it[iat];
             Atom* atom1 = &GlobalC::ucell.atoms[T1];
@@ -519,11 +569,13 @@ void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<
             // get iat1
             int iat1 = GlobalC::ucell.itia2iat(T1, I1);
             //
-            int irr = pv->nlocstart[iat];
+            int irr = pv.nlocstart[iat];
             const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0);
             double* ftvnl_dphi_iat;
-            if (isforce)
-                ftvnl_dphi_iat = &ftvnl_dphi(iat, 0);
+			if (isforce)
+			{
+				ftvnl_dphi_iat = &ftvnl_dphi(iat, 0);
+			}
 #ifdef _OPENMP
             // using local stack to avoid false sharing in multi-threaded case
             double ftvnl_dphi_temp[3] = {0.0, 0.0, 0.0};
@@ -544,7 +596,7 @@ void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<
                 double Ry = ra.info[iat][cb][1];
                 double Rz = ra.info[iat][cb][2];
                 // get BaseMatrix
-                if (pv->get_row_size(iat1) <= 0 || pv->get_col_size(iat2) <= 0)
+                if (pv.get_row_size(iat1) <= 0 || pv.get_col_size(iat2) <= 0)
                 {
                     continue;
                 }
@@ -554,9 +606,9 @@ void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<
                     tmp_matrix.push_back(DM->get_DMR_pointer(is+1)->find_matrix(iat1, iat2, Rx, Ry, Rz));
                 }
                 //hamilt::BaseMatrix<double>* tmp_matrix = DM->get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz);
-                for (int mu = 0; mu < pv->get_row_size(iat1); ++mu)
+                for (int mu = 0; mu < pv.get_row_size(iat1); ++mu)
                 {
-                    for (int nu = 0; nu < pv->get_col_size(iat2); ++nu)
+                    for (int nu = 0; nu < pv.get_col_size(iat2); ++nu)
                     {
                         // get value from DM
                         double dm2d1 = 0.0;
@@ -568,18 +620,18 @@ void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<
                         //
                         if (isforce)
                         {
-                            ftvnl_dphi_iat[0] += dm2d2 * this->UHM->LM->DHloc_fixedR_x[irr];
-                            ftvnl_dphi_iat[1] += dm2d2 * this->UHM->LM->DHloc_fixedR_y[irr];
-                            ftvnl_dphi_iat[2] += dm2d2 * this->UHM->LM->DHloc_fixedR_z[irr];
+                            ftvnl_dphi_iat[0] += dm2d2 * lm.DHloc_fixedR_x[irr];
+                            ftvnl_dphi_iat[1] += dm2d2 * lm.DHloc_fixedR_y[irr];
+                            ftvnl_dphi_iat[2] += dm2d2 * lm.DHloc_fixedR_z[irr];
                         }
                         if (isstress)
                         {
-                            local_stvnl_dphi(0, 0) -= dm2d1 * this->UHM->LM->stvnl11[irr];
-                            local_stvnl_dphi(0, 1) -= dm2d1 * this->UHM->LM->stvnl12[irr];
-                            local_stvnl_dphi(0, 2) -= dm2d1 * this->UHM->LM->stvnl13[irr];
-                            local_stvnl_dphi(1, 1) -= dm2d1 * this->UHM->LM->stvnl22[irr];
-                            local_stvnl_dphi(1, 2) -= dm2d1 * this->UHM->LM->stvnl23[irr];
-                            local_stvnl_dphi(2, 2) -= dm2d1 * this->UHM->LM->stvnl33[irr];
+                            local_stvnl_dphi(0, 0) -= dm2d1 * lm.stvnl11[irr];
+                            local_stvnl_dphi(0, 1) -= dm2d1 * lm.stvnl12[irr];
+                            local_stvnl_dphi(0, 2) -= dm2d1 * lm.stvnl13[irr];
+                            local_stvnl_dphi(1, 1) -= dm2d1 * lm.stvnl22[irr];
+                            local_stvnl_dphi(1, 2) -= dm2d1 * lm.stvnl23[irr];
+                            local_stvnl_dphi(2, 2) -= dm2d1 * lm.stvnl33[irr];
                         }
                         //}
                         ++local_total_irr;
@@ -612,10 +664,7 @@ void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<
         }
     }
 #endif
-    assert(total_irr == pv->nnr);
-
-    //	test(this->UHM->LM->DSloc_Rx);
-    //	test(dm2d[0],"dm2d");
+    assert(total_irr == pv.nnr);
 
     if (isstress)
     {
@@ -626,11 +675,17 @@ void Force_LCAO_k::cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<
     return;
 }
 
-void Force_LCAO_k::test(double* mmm, const std::string& name)
+void Force_LCAO_k::test(
+		Parallel_Orbitals &pv,
+		double* mmm, 
+		const std::string& name)
 {
-    const Parallel_Orbitals* pv = this->ParaV;
+    // mohan remove 'const' for pv, 2024-03-31
     if (GlobalV::NPROC != 1)
+    {
         return;
+    }
+
     std::cout << "test!" << std::endl;
 
     int irr = 0;
@@ -638,7 +693,9 @@ void Force_LCAO_k::test(double* mmm, const std::string& name)
 
     GlobalV::ofs_running << " Calculate the test in Force_LCAO_k" << std::endl;
     Record_adj RA;
-    RA.for_2d(*this->UHM->LM->ParaV, GlobalV::GAMMA_ONLY_LOCAL);
+    
+    // mohan update 2024-03-31
+    RA.for_2d(pv, GlobalV::GAMMA_ONLY_LOCAL);
 
     double* test;
     test = new double[GlobalV::NLOCAL * GlobalV::NLOCAL];
@@ -664,8 +721,7 @@ void Force_LCAO_k::test(double* mmm, const std::string& name)
                     for (int kk = 0; kk < atom2->nw; kk++)
                     {
                         const int iw2_all = start2 + kk;
-                        assert(irr < pv->nnr);
-                        // test[iw1_all*GlobalV::NLOCAL+iw2_all] += this->UHM->LM->DHloc_fixedR_x[irr];
+                        assert(irr < pv.nnr);
                         test[iw1_all * GlobalV::NLOCAL + iw2_all] += mmm[irr];
                         ++irr;
                     }
@@ -699,13 +755,13 @@ typedef std::tuple<int, int, int, int> key_tuple;
 // must consider three-center H matrix.
 void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<double>, double>* DM,
                                     const bool isforce,
-                                    const bool isstress,
-                                    ModuleBase::matrix& fvnl_dbeta,
+									const bool isstress,
+									const Parallel_Orbitals &pv,
+									ModuleBase::matrix& fvnl_dbeta,
                                     ModuleBase::matrix& svnl_dbeta)
 {
     ModuleBase::TITLE("Force_LCAO_k", "cal_fvnl_dbeta_k_new");
     ModuleBase::timer::tick("Force_LCAO_k", "cal_fvnl_dbeta_k_new");
-    const Parallel_Orbitals* pv = this->ParaV;
 
     // Data structure for storing <psi|beta>, for a detailed description
     // check out the same data structure in build_Nonlocal_mu_new
@@ -757,10 +813,12 @@ void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<
             for (int iw1 = 0; iw1 < nw1_tot; ++iw1)
             {
                 const int iw1_all = start1 + iw1;
-                const int iw1_local = pv->global2local_row(iw1_all);
-                const int iw2_local = pv->global2local_col(iw1_all);
-                if (iw1_local < 0 && iw2_local < 0)
-                    continue;
+                const int iw1_local = pv.global2local_row(iw1_all);
+                const int iw2_local = pv.global2local_col(iw1_all);
+				if (iw1_local < 0 && iw2_local < 0)
+				{
+					continue;
+				}
                 const int iw1_0 = iw1 / GlobalV::NPOL;
                 std::vector<std::vector<double>> nlm;
 #ifdef USE_NEW_TWO_CENTER
@@ -846,7 +904,7 @@ void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<
                 AdjacentAtomInfo adjs;
                 GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1, &adjs);
                 const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0);
-                int nnr = this->ParaV->nlocstart[iat1];
+                int nnr = pv.nlocstart[iat1];
 
                 /*
                     !!!!!!!!!!!!!!!!
@@ -894,8 +952,10 @@ void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<
                         for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0)
                         {
                             const int T0 = adjs.ntype[ad0];
-                            if (GlobalC::ucell.infoNL.nproj[T0] == 0)
-                                continue;
+							if (GlobalC::ucell.infoNL.nproj[T0] == 0)
+							{
+								continue;
+							}
                             const int I0 = adjs.natom[ad0];
                             // const int iat0 = GlobalC::ucell.itia2iat(T0, I0);
                             // const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0);
@@ -920,7 +980,7 @@ void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<
                     if (is_adj)
                     {
                         // basematrix and its data pointer
-                        if (pv->get_row_size(iat1) <= 0 || pv->get_col_size(iat2) <= 0)
+                        if (pv.get_row_size(iat1) <= 0 || pv.get_col_size(iat2) <= 0)
                         {
                             continue;
                         }
@@ -984,17 +1044,21 @@ void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<
                             {
                                 const int j0 = j / GlobalV::NPOL; // added by zhengdy-soc
                                 const int iw1_all = start1 + j;
-                                const int mu = pv->global2local_row(iw1_all);
-                                if (mu < 0)
-                                    continue;
+                                const int mu = pv.global2local_row(iw1_all);
+								if (mu < 0)
+								{
+									continue;
+								}
 
                                 for (int k = 0; k < atom2->nw * GlobalV::NPOL; k++)
                                 {
                                     const int k0 = k / GlobalV::NPOL;
                                     const int iw2_all = start2 + k;
-                                    const int nu = pv->global2local_col(iw2_all);
-                                    if (nu < 0)
-                                        continue;
+                                    const int nu = pv.global2local_col(iw2_all);
+									if (nu < 0)
+									{
+										continue;
+									}
 
                                     // const Atom* atom0 = &GlobalC::ucell.atoms[T0];
                                     double nlm[3] = {0, 0, 0};
@@ -1093,9 +1157,11 @@ void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<
                         {
                             const int j0 = j / GlobalV::NPOL; // added by zhengdy-soc
                             const int iw1_all = start1 + j;
-                            const int mu = pv->global2local_row(iw1_all);
-                            if (mu < 0)
-                                continue;
+                            const int mu = pv.global2local_row(iw1_all);
+							if (mu < 0)
+							{
+								continue;
+							}
 
                             // fix a serious bug: atom2[T2] -> atom2
                             // mohan 2010-12-20
@@ -1103,10 +1169,12 @@ void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<
                             {
                                 const int k0 = k / GlobalV::NPOL;
                                 const int iw2_all = start2 + k;
-                                const int nu = pv->global2local_col(iw2_all);
-                                if (nu < 0)
-                                    continue;
-                                total_nnr++;
+                                const int nu = pv.global2local_col(iw2_all);
+								if (nu < 0)
+								{
+									continue;
+								}
+								total_nnr++;
                                 nnr++;
                             }
                         }
@@ -1161,7 +1229,7 @@ void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<
     }
 #endif
 
-    assert(total_nnr == pv->nnr);
+    assert(total_nnr == pv.nnr);
 
     if (isstress)
     {
@@ -1174,11 +1242,12 @@ void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<
 
 // calculate the force due to < phi | Vlocal | dphi >
 void Force_LCAO_k::cal_fvl_dphi_k(const bool isforce,
-                                  const bool isstress,
-                                  const elecstate::Potential* pot_in,
-                                  ModuleBase::matrix& fvl_dphi,
-                                  ModuleBase::matrix& svl_dphi,
-                                  double** DM_R)
+		const bool isstress,
+		LCAO_Matrix &lm,
+		const elecstate::Potential* pot_in,
+		ModuleBase::matrix& fvl_dphi,
+		ModuleBase::matrix& svl_dphi,
+		double** DM_R)
 {
     ModuleBase::TITLE("Force_LCAO_k", "cal_fvl_dphi_k");
     ModuleBase::timer::tick("Force_LCAO_k", "cal_fvl_dphi_k");
@@ -1188,9 +1257,10 @@ void Force_LCAO_k::cal_fvl_dphi_k(const bool isforce,
         ModuleBase::timer::tick("Force_LCAO_k", "cal_fvl_dphi_k");
         return;
     }
-    assert(this->UHM->LM->DHloc_fixedR_x != NULL);
-    assert(this->UHM->LM->DHloc_fixedR_y != NULL);
-    assert(this->UHM->LM->DHloc_fixedR_z != NULL);
+
+    assert(lm.DHloc_fixedR_x != NULL);
+    assert(lm.DHloc_fixedR_y != NULL);
+    assert(lm.DHloc_fixedR_z != NULL);
 
     int istep = 1;
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h
index 0deb189aa56..2d9b3bc2c30 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.h
@@ -12,16 +12,15 @@
 
 class Force_LCAO_k : public Force_LCAO_gamma
 {
-public:
+    public:
     template <typename T>
     friend class Force_Stress_LCAO;
 
     Force_LCAO_k();
     ~Force_LCAO_k();
 
-  private:
+    private:
     LCAO_Hamilt* UHM;
-    const Parallel_Orbitals* ParaV;
 
     // orthonormal force + contribution from T and VNL
     void ftable_k(const bool isforce,
@@ -43,53 +42,65 @@ class Force_LCAO_k : public Force_LCAO_gamma
 #else
                   ModuleBase::matrix& svl_dphi,
 #endif
-                  LCAO_Hamilt& uhm,
-                  const K_Vectors& kv);
+				  LCAO_Hamilt &uhm,
+				  Parallel_Orbitals &pv,
+				  LCAO_Matrix &lm,
+                  const K_Vectors &kv);
 
     // get the ds, dt, dvnl.
     void allocate_k(const Parallel_Orbitals& pv,
+                    LCAO_Matrix &lm,
                     const int& nks,
                     const std::vector<ModuleBase::Vector3<double>>& kvec_d);
 
-    void finish_k(void);
+    void finish_k(LCAO_Matrix &lm);
 
     // calculate the force due to < dphi | beta > < beta | phi >
-    void cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<double>, double>* DM,
-                          const bool isforce,
-                          const bool isstress,
-                          Record_adj& ra,
-                          ModuleBase::matrix& ftvnl_dphi,
-                          ModuleBase::matrix& stvnl_dphi);
+	void cal_ftvnl_dphi_k(const elecstate::DensityMatrix<std::complex<double>, double>* DM,
+			const Parallel_Orbitals &pv,
+			LCAO_Matrix &lm,
+			const bool isforce,
+			const bool isstress,
+			Record_adj& ra,
+			ModuleBase::matrix& ftvnl_dphi,
+			ModuleBase::matrix& stvnl_dphi);
 
     // calculate the overlap force
     void cal_foverlap_k(const bool isforce,
                         const bool isstress,
-                        Record_adj& ra,
-                        const psi::Psi<std::complex<double>>* psi,
-                        Local_Orbital_Charge& loc,
-                        const elecstate::DensityMatrix<std::complex<double>, double>* DM,
-                        ModuleBase::matrix& foverlap,
-                        ModuleBase::matrix& soverlap,
-                        const elecstate::ElecState* pelec,
-                        const int& nks,
-                        const K_Vectors& kv);
+                        Record_adj &ra,
+                        const psi::Psi<std::complex<double>> *psi,
+						Local_Orbital_Charge &loc,
+						Parallel_Orbitals &pv,
+						LCAO_Matrix &lm,
+						const elecstate::DensityMatrix<std::complex<double>, double> *DM,
+                        ModuleBase::matrix &foverlap,
+                        ModuleBase::matrix &soverlap,
+                        const elecstate::ElecState *pelec,
+                        const int &nks,
+                        const K_Vectors &kv);
 
     // calculate the force due to < phi | Vlocal | dphi >
-    void cal_fvl_dphi_k(const bool isforce,
-                        const bool isstress,
-                        const elecstate::Potential* pot_in,
-                        ModuleBase::matrix& fvl_dphi,
-                        ModuleBase::matrix& svl_dphi,
-                        double** DM_R);
+	void cal_fvl_dphi_k(const bool isforce,
+		  	const bool isstress,
+			LCAO_Matrix &lm,
+			const elecstate::Potential* pot_in,
+			ModuleBase::matrix& fvl_dphi,
+			ModuleBase::matrix& svl_dphi,
+			double** DM_R);
 
     // new method to calculate the force due to < phi | dbeta > < beta | phi > , developed by wenfei-li
     void cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<double>, double>* DM,
                           const bool isforce,
-                          const bool isstress, 
-                          ModuleBase::matrix& fvnl_dbeta,
+						  const bool isstress, 
+						  const Parallel_Orbitals &pv,
+						  ModuleBase::matrix& fvnl_dbeta,
                           ModuleBase::matrix& svnl_dbeta);
 
-    void test(double* mm, const std::string& name);
+	void test(
+			Parallel_Orbitals &pv,
+			double* mm, 
+			const std::string& name);
 };
 
 #endif
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp
index 4030ff5eef3..419f1726e71 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp
@@ -31,7 +31,10 @@ LCAO_Hamilt::~LCAO_Hamilt()
 //--------------------------------------------
 // prepare grid network for Gint(grid integral)
 //--------------------------------------------
-void LCAO_Hamilt::grid_prepare(const Grid_Technique& gt, const ModulePW::PW_Basis& rhopw, const ModulePW::PW_Basis_Big& bigpw)
+void LCAO_Hamilt::grid_prepare(
+		const Grid_Technique& gt, 
+		const ModulePW::PW_Basis& rhopw, 
+		const ModulePW::PW_Basis_Big& bigpw)
 {
     ModuleBase::TITLE("LCAO_Hamilt","grid_prepare");
     ModuleBase::timer::tick("LCAO_Hamilt","grid_prepare");
@@ -45,7 +48,7 @@ void LCAO_Hamilt::grid_prepare(const Grid_Technique& gt, const ModulePW::PW_Basi
     }
     else // multiple k-points
     {
-        // calculate the grid integration of 'Vl' matrix for l-points algorithms.
+        // cal the grid integration of 'Vl' matrix for l-points algorithms.
         this->GK.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);
@@ -55,7 +58,7 @@ void LCAO_Hamilt::grid_prepare(const Grid_Technique& gt, const ModulePW::PW_Basi
     return;
 }
 
-void LCAO_Hamilt::set_R_range_sparse()
+void LCAO_Hamilt::set_R_range_sparse(void)
 {
     int R_minX = int(GlobalC::GridD.getD_minX());
     int R_minY = int(GlobalC::GridD.getD_minY());
@@ -80,15 +83,15 @@ void LCAO_Hamilt::set_R_range_sparse()
     return;
 }
 
-void LCAO_Hamilt::calculate_STN_R_sparse_for_S(const double &sparse_threshold)
+void LCAO_Hamilt::cal_STN_R_sparse_for_S(const double &sparse_threshold)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculate_STN_R_sparse_for_S");
+    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 temp_value_double;
+    double temp_value_double = 0.0;
     std::complex<double> temp_value_complex;
 
     for(int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1)
@@ -188,9 +191,13 @@ void LCAO_Hamilt::calculate_STN_R_sparse_for_S(const double &sparse_threshold)
 }
 
 #include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
-void LCAO_Hamilt::calculate_HContainer_sparse_d(const int &current_spin, const double &sparse_threshold, const hamilt::HContainer<double>& hR, std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>>& target)
+void LCAO_Hamilt::cal_HContainer_sparse_d(
+		const int &current_spin, 
+		const double &sparse_threshold, 
+		const hamilt::HContainer<double>& hR, 
+		std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>>& target)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculate_HContainer_sparse_d");
+    ModuleBase::TITLE("LCAO_Hamilt","cal_HContainer_sparse_d");
 
     const Parallel_Orbitals* paraV = this->LM->ParaV;
     auto row_indexes = paraV->get_indexes_row();
@@ -227,9 +234,14 @@ void LCAO_Hamilt::calculate_HContainer_sparse_d(const int &current_spin, const d
     return;
 }
 
-void LCAO_Hamilt::calculate_HContainer_sparse_cd(const int &current_spin, const double &sparse_threshold, 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 LCAO_Hamilt::cal_HContainer_sparse_cd(
+		const int &current_spin, 
+		const double &sparse_threshold, 
+		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)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculate_HContainer_sparse_cd");
+    ModuleBase::TITLE("LCAO_Hamilt","cal_HContainer_sparse_cd");
 
     const Parallel_Orbitals* paraV = this->LM->ParaV;
     auto row_indexes = paraV->get_indexes_row();
@@ -266,36 +278,58 @@ void LCAO_Hamilt::calculate_HContainer_sparse_cd(const int &current_spin, const
     return;
 }
 
-void LCAO_Hamilt::calculate_HSR_sparse(const int &current_spin, const double &sparse_threshold, const int (&nmp)[3], hamilt::Hamilt<std::complex<double>>* p_ham)
+void LCAO_Hamilt::cal_HSR_sparse(
+		const int &current_spin, 
+		const double &sparse_threshold, 
+		const int (&nmp)[3], 
+		hamilt::Hamilt<std::complex<double>>* p_ham)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculate_HSR_sparse");
+    ModuleBase::TITLE("LCAO_Hamilt","cal_HSR_sparse");
 
     set_R_range_sparse();
 
-    //calculate_STN_R_sparse(current_spin, sparse_threshold);
+    //cal_STN_R_sparse(current_spin, sparse_threshold);
     if(GlobalV::NSPIN!=4)
     {
-        hamilt::HamiltLCAO<std::complex<double>, double>* p_ham_lcao = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, double>*>(p_ham);
-        this->calculate_HContainer_sparse_d(current_spin, sparse_threshold, *(p_ham_lcao->getHR()), this->LM->HR_sparse[current_spin]);
-        this->calculate_HContainer_sparse_d(current_spin, sparse_threshold, *(p_ham_lcao->getSR()), this->LM->SR_sparse);
+        hamilt::HamiltLCAO<std::complex<double>, double>* p_ham_lcao = 
+        dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, double>*>(p_ham);
+
+		this->cal_HContainer_sparse_d(current_spin, 
+				sparse_threshold, 
+				*(p_ham_lcao->getHR()), 
+				this->LM->HR_sparse[current_spin]);
+
+		this->cal_HContainer_sparse_d(current_spin, 
+				sparse_threshold, 
+				*(p_ham_lcao->getSR()), 
+				this->LM->SR_sparse);
     }
     else
     {
-        hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>* p_ham_lcao = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>*>(p_ham);
-        this->calculate_HContainer_sparse_cd(current_spin, sparse_threshold, *(p_ham_lcao->getHR()), this->LM->HR_soc_sparse);
-        this->calculate_HContainer_sparse_cd(current_spin, sparse_threshold, *(p_ham_lcao->getSR()), this->LM->SR_soc_sparse);
+        hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>* p_ham_lcao = 
+        dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>*>(p_ham);
+
+        this->cal_HContainer_sparse_cd(current_spin, 
+        sparse_threshold, 
+        *(p_ham_lcao->getHR()), 
+        this->LM->HR_soc_sparse);
+
+        this->cal_HContainer_sparse_cd(current_spin, 
+        sparse_threshold, 
+        *(p_ham_lcao->getSR()), 
+        this->LM->SR_soc_sparse);
     }
 
-    // only old DFT+U method need to calculate extra contribution to HR
+    // only old DFT+U method need to cal extra contribution to HR
     if (GlobalV::dft_plus_u == 2)
     {
         if (GlobalV::NSPIN != 4)
         {
-            calculat_HR_dftu_sparse(current_spin, sparse_threshold);
+            cal_HR_dftu_sparse(current_spin, sparse_threshold);
         }
         else
         {
-            calculat_HR_dftu_soc_sparse(current_spin, sparse_threshold);
+            cal_HR_dftu_soc_sparse(current_spin, sparse_threshold);
         }
     }
 
@@ -303,20 +337,24 @@ void LCAO_Hamilt::calculate_HSR_sparse(const int &current_spin, const double &sp
 #ifdef __MPI
     if( GlobalC::exx_info.info_global.cal_exx )
     {
-        if(GlobalC::exx_info.info_ri.real_number)
-            this->calculate_HR_exx_sparse(current_spin, sparse_threshold, nmp, *this->LM->Hexxd);
-        else
-            this->calculate_HR_exx_sparse(current_spin, sparse_threshold, nmp, *this->LM->Hexxc);
-    }
+		if(GlobalC::exx_info.info_ri.real_number)
+		{
+			this->cal_HR_exx_sparse(current_spin, sparse_threshold, nmp, *this->LM->Hexxd);
+		}
+		else
+		{
+			this->cal_HR_exx_sparse(current_spin, sparse_threshold, nmp, *this->LM->Hexxc);
+		}
+	}
 #endif // __MPI
 #endif // __EXX
 
     clear_zero_elements(current_spin, sparse_threshold);
 }
 
-void LCAO_Hamilt::calculate_dH_sparse(const int &current_spin, const double &sparse_threshold)
+void LCAO_Hamilt::cal_dH_sparse(const int &current_spin, const double &sparse_threshold)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculate_dH_sparse");
+    ModuleBase::TITLE("LCAO_Hamilt","cal_dH_sparse");
 
     set_R_range_sparse();
 
@@ -328,8 +366,8 @@ void LCAO_Hamilt::calculate_dH_sparse(const int &current_spin, const double &spa
     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);
-    // calculate dT=<phi|kin|dphi> in LCAO
-    // calculate T + VNL(P1) in LCAO basis
+    // cal dT=<phi|kin|dphi> in LCAO
+    // cal T + VNL(P1) in LCAO basis
     if(GlobalV::CAL_STRESS)
 	{
         GlobalV::CAL_STRESS = false;
@@ -342,7 +380,7 @@ void LCAO_Hamilt::calculate_dH_sparse(const int &current_spin, const double &spa
     }
     this->genH.build_Nonlocal_mu_new (this->LM->Hloc_fixed.data(), true);
     
-    calculate_dSTN_R_sparse(current_spin, sparse_threshold);
+    cal_dSTN_R_sparse(current_spin, sparse_threshold);
 
     delete[] this->LM->DHloc_fixedR_x;
     delete[] this->LM->DHloc_fixedR_y;
@@ -351,9 +389,9 @@ void LCAO_Hamilt::calculate_dH_sparse(const int &current_spin, const double &spa
     GK.cal_dvlocal_R_sparseMatrix(current_spin, sparse_threshold, this->LM);
 }
 
-void LCAO_Hamilt::calculate_STN_R_sparse_for_T(const double &sparse_threshold)
+void LCAO_Hamilt::cal_STN_R_sparse_for_T(const double &sparse_threshold)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculate_STN_R_sparse_for_T");
+    ModuleBase::TITLE("LCAO_Hamilt","cal_STN_R_sparse_for_T");
 
     int index = 0;
     ModuleBase::Vector3<double> dtau, tau1, tau2;
@@ -414,7 +452,10 @@ void LCAO_Hamilt::calculate_STN_R_sparse_for_T(const double &sparse_threshold)
                 {
                     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);
+					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++)
                     {
@@ -458,39 +499,43 @@ void LCAO_Hamilt::calculate_STN_R_sparse_for_T(const double &sparse_threshold)
     return;
 }
 
-void LCAO_Hamilt::calculate_SR_sparse(const double &sparse_threshold, hamilt::Hamilt<std::complex<double>>* p_ham)
+void LCAO_Hamilt::cal_SR_sparse(const double &sparse_threshold, hamilt::Hamilt<std::complex<double>>* p_ham)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculate_SR_sparse");
+    ModuleBase::TITLE("LCAO_Hamilt","cal_SR_sparse");
     set_R_range_sparse();
-    //calculate_STN_R_sparse(current_spin, sparse_threshold);
+    //cal_STN_R_sparse(current_spin, sparse_threshold);
     if(GlobalV::NSPIN!=4)
     {
-        hamilt::HamiltLCAO<std::complex<double>, double>* p_ham_lcao = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, double>*>(p_ham);
-        this->calculate_HContainer_sparse_d(0, sparse_threshold, *(p_ham_lcao->getSR()), this->LM->SR_sparse);
+        hamilt::HamiltLCAO<std::complex<double>, double>* p_ham_lcao 
+        = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, double>*>(p_ham);
+        this->cal_HContainer_sparse_d(0, sparse_threshold, *(p_ham_lcao->getSR()), this->LM->SR_sparse);
     }
     else
     {
-        hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>* p_ham_lcao = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>*>(p_ham);
-        this->calculate_HContainer_sparse_cd(0, sparse_threshold, *(p_ham_lcao->getSR()), this->LM->SR_soc_sparse);
+        hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>* p_ham_lcao 
+        = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>*>(p_ham);
+        this->cal_HContainer_sparse_cd(0, sparse_threshold, *(p_ham_lcao->getSR()), this->LM->SR_soc_sparse);
     }
 }
 
-void LCAO_Hamilt::calculate_TR_sparse(const double &sparse_threshold)
+void LCAO_Hamilt::cal_TR_sparse(const double &sparse_threshold)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculate_TR_sparse");
+    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());
-    set_R_range_sparse();
-    calculate_STN_R_sparse_for_T(sparse_threshold);
+    this->set_R_range_sparse();
+    this->cal_STN_R_sparse_for_T(sparse_threshold);
+
+    return;
 }
 
-void LCAO_Hamilt::calculat_HR_dftu_sparse(const int &current_spin, const double &sparse_threshold)
+void LCAO_Hamilt::cal_HR_dftu_sparse(const int &current_spin, const double &sparse_threshold)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculat_HR_dftu_sparse");
-    ModuleBase::timer::tick("LCAO_Hamilt","calculat_HR_dftu_sparse");
+    ModuleBase::TITLE("LCAO_Hamilt","cal_HR_dftu_sparse");
+    ModuleBase::timer::tick("LCAO_Hamilt","cal_HR_dftu_sparse");
 
     int total_R_num = this->LM->all_R_coor.size();
     int *nonzero_num = new int[total_R_num];
@@ -514,9 +559,9 @@ void LCAO_Hamilt::calculat_HR_dftu_sparse(const int &current_spin, const double
     double *HR_tmp = new double[this->LM->ParaV->nloc];
     double *SR_tmp = new double[this->LM->ParaV->nloc];
 
-    int ir;
-    int ic;
-    int iic;
+    int ir=0;
+    int ic=0;
+    int iic=0;
     auto &temp_HR_sparse = this->LM->HR_sparse[current_spin];
 
     count = 0;
@@ -596,14 +641,14 @@ void LCAO_Hamilt::calculat_HR_dftu_sparse(const int &current_spin, const double
     HR_tmp = nullptr;
     SR_tmp = nullptr;
 
-    ModuleBase::timer::tick("LCAO_Hamilt","calculat_HR_dftu_sparse");
+    ModuleBase::timer::tick("LCAO_Hamilt","cal_HR_dftu_sparse");
 
 }
 
-void LCAO_Hamilt::calculat_HR_dftu_soc_sparse(const int &current_spin, const double &sparse_threshold)
+void LCAO_Hamilt::cal_HR_dftu_soc_sparse(const int &current_spin, const double &sparse_threshold)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculat_HR_dftu_soc_sparse");
-    ModuleBase::timer::tick("LCAO_Hamilt","calculat_HR_dftu_soc_sparse");
+    ModuleBase::TITLE("LCAO_Hamilt","cal_HR_dftu_soc_sparse");
+    ModuleBase::timer::tick("LCAO_Hamilt","cal_HR_dftu_soc_sparse");
 
     int total_R_num = this->LM->all_R_coor.size();
     int *nonzero_num = new int[total_R_num];
@@ -627,9 +672,9 @@ void LCAO_Hamilt::calculat_HR_dftu_soc_sparse(const int &current_spin, const dou
     std::complex<double> *HR_soc_tmp = new std::complex<double>[this->LM->ParaV->nloc];
     std::complex<double> *SR_soc_tmp = new std::complex<double>[this->LM->ParaV->nloc];
 
-    int ir;
-    int ic;
-    int iic;
+    int ir=0;
+    int ic=0;
+    int iic=0;
 
     count = 0;
     for (auto &R_coor : this->LM->all_R_coor)
@@ -803,9 +848,9 @@ void LCAO_Hamilt::clear_zero_elements(const int &current_spin, const double &spa
     }
 }
 
-void LCAO_Hamilt::calculate_dSTN_R_sparse(const int &current_spin, const double &sparse_threshold)
+void LCAO_Hamilt::cal_dSTN_R_sparse(const int &current_spin, const double &sparse_threshold)
 {
-    ModuleBase::TITLE("LCAO_Hamilt","calculate_dSTN_R_sparse");
+    ModuleBase::TITLE("LCAO_Hamilt","cal_dSTN_R_sparse");
 
     int index = 0;
     ModuleBase::Vector3<double> dtau, tau1, tau2;
@@ -902,7 +947,7 @@ void LCAO_Hamilt::calculate_dSTN_R_sparse(const int &current_spin, const double
                             }
                             else
                             {
-                                ModuleBase::WARNING_QUIT("calculate_dSTN_R_sparse","soc not supported!");
+                                ModuleBase::WARNING_QUIT("cal_dSTN_R_sparse","soc not supported!");
                             }
                             ++index;
                         }
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h
index 90d23c4be30..1422752c15d 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h
@@ -26,51 +26,51 @@ class LCAO_Hamilt
     // jingan add 2021-6-4
     void set_R_range_sparse();
 
-    void calculate_HContainer_sparse_d(const int &current_spin, 
+    void cal_HContainer_sparse_d(const int &current_spin, 
         const double &sparse_threshold, 
         const hamilt::HContainer<double>& hR, 
         std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>>& target);
 
-    void calculate_HContainer_sparse_cd(const int &current_spin, 
+    void cal_HContainer_sparse_cd(const int &current_spin, 
         const double &sparse_threshold, 
         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 calculate_dSTN_R_sparse(const int &current_spin, const double &sparse_threshold);
+    void cal_dSTN_R_sparse(const int &current_spin, const double &sparse_threshold);
 
-    void calculate_STN_R_sparse_for_S(const double &sparse_threshold);
+    void cal_STN_R_sparse_for_S(const double &sparse_threshold);
 
-    void calculate_STN_R_sparse_for_T(const double &sparse_threshold);
+    void cal_STN_R_sparse_for_T(const double &sparse_threshold);
 
-    void calculat_HR_dftu_sparse(const int &current_spin, const double &sparse_threshold);
+    void cal_HR_dftu_sparse(const int &current_spin, const double &sparse_threshold);
 
-    void calculat_HR_dftu_soc_sparse(const int &current_spin, const double &sparse_threshold);
+    void cal_HR_dftu_soc_sparse(const int &current_spin, const double &sparse_threshold);
 
 #ifdef __EXX
-    template<typename Tdata> void calculate_HR_exx_sparse(
+    template<typename Tdata> void cal_HR_exx_sparse(
             const int &current_spin,
             const double &sparse_threshold,
             const int (&nmp)[3],
             const std::vector< std::map <int, std::map < std::pair<int, std::array<int,3>>, RI::Tensor<Tdata> > >>& Hexxs);
 #endif
 
-	void calculate_HSR_sparse(
+	void cal_HSR_sparse(
 			const int &current_spin, 
 			const double &sparse_threshold, 
 			const int (&nmp)[3], 
 			hamilt::Hamilt<std::complex<double>>* p_ham);
 
-    void calculate_SR_sparse(const double &sparse_threshold, hamilt::Hamilt<std::complex<double>>* p_ham);
+    void cal_SR_sparse(const double &sparse_threshold, hamilt::Hamilt<std::complex<double>>* p_ham);
 
     void clear_zero_elements(const int &current_spin, const double &sparse_threshold);
 
     void destroy_all_HSR_sparse(void);
 
-    void calculate_TR_sparse(const double &sparse_threshold);
+    void cal_TR_sparse(const double &sparse_threshold);
 
     void destroy_TR_sparse(void);
 
-    void calculate_dH_sparse(const int &current_spin, const double &sparse_threshold);
+    void cal_dH_sparse(const int &current_spin, const double &sparse_threshold);
 
     void destroy_dH_R_sparse(void);
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp
index db8b32a13ab..57ed48c7e7f 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp
@@ -23,30 +23,35 @@
 
 // Peize Lin add 2022.09.13
 template<typename Tdata>
-void LCAO_Hamilt::calculate_HR_exx_sparse(
+void LCAO_Hamilt::cal_HR_exx_sparse(
 			const int &current_spin, 
 			const double &sparse_threshold,
 			const int (&nmp)[3],
 			const std::vector< std::map<int, std::map<std::pair<int,std::array<int,3>>, RI::Tensor<Tdata>>>>& Hexxs)
 {
-	ModuleBase::TITLE("LCAO_Hamilt","calculate_HR_exx_sparse");
-	ModuleBase::timer::tick("LCAO_Hamilt","calculate_HR_exx_sparse");
+	ModuleBase::TITLE("LCAO_Hamilt","cal_HR_exx_sparse");
+	ModuleBase::timer::tick("LCAO_Hamilt","cal_HR_exx_sparse");
 
 	const Tdata frac = GlobalC::exx_info.info_global.hybrid_alpha;
 
 	std::map<int,std::array<double,3>> atoms_pos;
 	for(int iat=0; iat<GlobalC::ucell.nat; ++iat)
-		atoms_pos[iat] = RI_Util::Vector3_to_array3( GlobalC::ucell.atoms[ GlobalC::ucell.iat2it[iat] ].tau[ GlobalC::ucell.iat2ia[iat] ] );
+	{
+		atoms_pos[iat] = RI_Util::Vector3_to_array3( 
+        GlobalC::ucell.atoms[ GlobalC::ucell.iat2it[iat] ].tau[ GlobalC::ucell.iat2ia[iat] ] );
+	}
 	const std::array<std::array<double,3>,3> latvec
 		= {RI_Util::Vector3_to_array3(GlobalC::ucell.a1),
 		   RI_Util::Vector3_to_array3(GlobalC::ucell.a2),
 		   RI_Util::Vector3_to_array3(GlobalC::ucell.a3)};
+
 	const std::array<int,3> Rs_period = {nmp[0], nmp[1], nmp[2]};
 
 	RI::Cell_Nearest<int,int,3,double,3> cell_nearest;
 	cell_nearest.init(atoms_pos, latvec, Rs_period);	
 
 	const std::vector<int> is_list = (GlobalV::NSPIN!=4) ? std::vector<int>{current_spin} : std::vector<int>{0,1,2,3};
+
 	for(const int is : is_list)
     {
 		int is0_b, is1_b;
@@ -102,7 +107,7 @@ void LCAO_Hamilt::calculate_HR_exx_sparse(
 		}
 	}
 
-	ModuleBase::timer::tick("LCAO_Hamilt","calculate_HR_exx_sparse");
+	ModuleBase::timer::tick("LCAO_Hamilt","cal_HR_exx_sparse");
 }
 
 #endif
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp
new file mode 100644
index 00000000000..6cfac40f2d7
--- /dev/null
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp
@@ -0,0 +1,489 @@
+
+// must consider three-center H matrix.
+void Force_LCAO_k::cal_fvnl_dbeta_k(const elecstate::DensityMatrix<std::complex<double>, double>* DM,
+                                    const bool isforce,
+									const bool isstress,
+									const Parallel_Orbitals &pv,
+									ModuleBase::matrix& fvnl_dbeta,
+                                    ModuleBase::matrix& svnl_dbeta)
+{
+    ModuleBase::TITLE("Force_LCAO_k", "cal_fvnl_dbeta_k_new");
+    ModuleBase::timer::tick("Force_LCAO_k", "cal_fvnl_dbeta_k_new");
+
+    // Data structure for storing <psi|beta>, for a detailed description
+    // check out the same data structure in build_Nonlocal_mu_new
+    std::vector<std::map<key_tuple, std::unordered_map<int, std::vector<std::vector<double>>>>> nlm_tot;
+
+    nlm_tot.resize(GlobalC::ucell.nat);
+
+#ifdef _OPENMP
+// use schedule(dynamic) for load balancing because adj_num is various
+#pragma omp parallel for schedule(dynamic)
+#endif
+    for (int iat = 0; iat < GlobalC::ucell.nat; iat++)
+    {
+
+        const int it = GlobalC::ucell.iat2it[iat];
+        const int ia = GlobalC::ucell.iat2ia[iat];
+
+        // Step 1 : generate <psi|beta>
+        // type of atom; distance; atomic basis; projectors
+
+        const double Rcut_Beta = GlobalC::ucell.infoNL.Beta[it].get_rcut_max();
+        const ModuleBase::Vector3<double> tau = GlobalC::ucell.atoms[it].tau[ia];
+        AdjacentAtomInfo adjs;
+        GlobalC::GridD.Find_atom(GlobalC::ucell, tau, it, ia, &adjs);
+
+        nlm_tot[iat].clear();
+
+        for (int ad = 0; ad < adjs.adj_num + 1; ++ad)
+        {
+            const int T1 = adjs.ntype[ad];
+            const int I1 = adjs.natom[ad];
+            const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0);
+            const double Rcut_AO1 = GlobalC::ORB.Phi[T1].getRcut();
+
+            const ModuleBase::Vector3<double>& tau1 = adjs.adjacent_tau[ad];
+            const Atom* atom1 = &GlobalC::ucell.atoms[T1];
+            const int nw1_tot = atom1->nw * GlobalV::NPOL;
+
+            const ModuleBase::Vector3<double> dtau = tau1 - tau;
+            const double dist1 = dtau.norm2() * pow(GlobalC::ucell.lat0, 2);
+            if (dist1 > pow(Rcut_Beta + Rcut_AO1, 2))
+            {
+                continue;
+            }
+
+            std::unordered_map<int, std::vector<std::vector<double>>> nlm_cur;
+            nlm_cur.clear();
+
+            for (int iw1 = 0; iw1 < nw1_tot; ++iw1)
+            {
+                const int iw1_all = start1 + iw1;
+                const int iw1_local = pv.global2local_row(iw1_all);
+                const int iw2_local = pv.global2local_col(iw1_all);
+				if (iw1_local < 0 && iw2_local < 0)
+				{
+					continue;
+				}
+                const int iw1_0 = iw1 / GlobalV::NPOL;
+                std::vector<std::vector<double>> nlm;
+#ifdef USE_NEW_TWO_CENTER
+                //=================================================================
+                //          new two-center integral (temporary)
+                //=================================================================
+                int L1 = atom1->iw2l[ iw1_0 ];
+                int N1 = atom1->iw2n[ iw1_0 ];
+                int m1 = atom1->iw2m[ iw1_0 ];
+
+                // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l)
+                int M1 = (m1 % 2 == 0) ? -m1/2 : (m1+1)/2;
+
+                ModuleBase::Vector3<double> dtau = tau - tau1;
+                GlobalC::UOT.two_center_bundle->overlap_orb_beta->snap(
+                        T1, L1, N1, M1, it, dtau * GlobalC::ucell.lat0, true, nlm);
+#else
+                GlobalC::UOT.snap_psibeta_half(GlobalC::ORB,
+                                               GlobalC::ucell.infoNL,
+                                               nlm,
+                                               tau1,
+                                               T1,
+                                               atom1->iw2l[iw1_0], // L1
+                                               atom1->iw2m[iw1_0], // m1
+                                               atom1->iw2n[iw1_0], // N1
+                                               tau,
+                                               it,
+                                               1); // R0,T0
+#endif
+                nlm_cur.insert({iw1_all, nlm});
+            } // end iw
+            const int iat1 = GlobalC::ucell.itia2iat(T1, I1);
+            const int rx1 = adjs.box[ad].x;
+            const int ry1 = adjs.box[ad].y;
+            const int rz1 = adjs.box[ad].z;
+            key_tuple key_1(iat1, rx1, ry1, rz1);
+            nlm_tot[iat][key_1] = nlm_cur;
+        } // end ad
+    }
+
+    //=======================================================
+    // Step2:
+    // calculate sum_(L0,M0) beta<psi_i|beta><beta|psi_j>
+    // and accumulate the value to Hloc_fixedR(i,j)
+    //=======================================================
+    int total_nnr = 0;
+#ifdef _OPENMP
+#pragma omp parallel reduction(+ : total_nnr)
+    {
+		ModuleBase::matrix local_svnl_dbeta(3, 3);
+		const int num_threads = omp_get_num_threads();
+#else
+		ModuleBase::matrix& local_svnl_dbeta = svnl_dbeta;
+#endif
+
+        ModuleBase::Vector3<double> tau1;
+        ModuleBase::Vector3<double> tau2;
+        ModuleBase::Vector3<double> dtau;
+        ModuleBase::Vector3<double> tau0;
+        ModuleBase::Vector3<double> dtau1;
+        ModuleBase::Vector3<double> dtau2;
+
+        double rcut;
+        double distance;
+
+        double rcut1;
+        double rcut2;
+        double distance1;
+        double distance2;
+
+#ifdef _OPENMP
+// use schedule(dynamic) for load balancing because adj_num is various
+#pragma omp for schedule(dynamic)
+#endif
+        for (int iat1 = 0; iat1 < GlobalC::ucell.nat; iat1++)
+        {
+            const int T1 = GlobalC::ucell.iat2it[iat1];
+            const Atom* atom1 = &GlobalC::ucell.atoms[T1];
+
+            {
+                const int I1 = GlobalC::ucell.iat2ia[iat1];
+                tau1 = atom1->tau[I1];
+                AdjacentAtomInfo adjs;
+                GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1, &adjs);
+                const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0);
+                int nnr = pv.nlocstart[iat1];
+
+                /*
+                    !!!!!!!!!!!!!!!!
+                    This optimization is also improving the performance of single thread.
+                    Making memory access more linearly in the core loop
+                */
+                bool iat_recorded = false;
+                bool force_updated = false;
+                // record iat of adjs
+                std::vector<int> adj_iat;
+                // record fvnl_dbeta diff of adjs
+                std::vector<double> adj_fvnl_dbeta;
+                if (isforce)
+                {
+                    adj_iat.resize(adjs.adj_num + 1);
+                    adj_fvnl_dbeta.resize((adjs.adj_num + 1) * 3, 0.0);
+                }
+
+                for (int ad2 = 0; ad2 < adjs.adj_num + 1; ++ad2)
+                {
+                    const int T2 = adjs.ntype[ad2];
+                    const Atom* atom2 = &GlobalC::ucell.atoms[T2];
+
+                    const int I2 = adjs.natom[ad2];
+                    const int iat2 = GlobalC::ucell.itia2iat(T2, I2);
+                    const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0);
+                    tau2 = adjs.adjacent_tau[ad2];
+
+                    const int rx2 = adjs.box[ad2].x;
+                    const int ry2 = adjs.box[ad2].y;
+                    const int rz2 = adjs.box[ad2].z;
+
+                    dtau = tau2 - tau1;
+                    distance = dtau.norm2() * pow(GlobalC::ucell.lat0, 2);
+                    // this rcut is in order to make nnr consistent
+                    // with other matrix.
+                    rcut = pow(GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ORB.Phi[T2].getRcut(), 2);
+
+                    // check if this a adjacent atoms.
+                    bool is_adj = false;
+                    if (distance < rcut)
+                        is_adj = true;
+                    else if (distance >= rcut)
+                    {
+                        for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0)
+                        {
+                            const int T0 = adjs.ntype[ad0];
+							if (GlobalC::ucell.infoNL.nproj[T0] == 0)
+							{
+								continue;
+							}
+                            const int I0 = adjs.natom[ad0];
+                            // const int iat0 = GlobalC::ucell.itia2iat(T0, I0);
+                            // const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0);
+
+                            tau0 = adjs.adjacent_tau[ad0];
+                            dtau1 = tau0 - tau1;
+                            distance1 = dtau1.norm() * GlobalC::ucell.lat0;
+                            rcut1 = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
+
+                            dtau2 = tau0 - tau2;
+                            distance2 = dtau2.norm() * GlobalC::ucell.lat0;
+                            rcut2 = GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
+
+                            if (distance1 < rcut1 && distance2 < rcut2)
+                            {
+                                is_adj = true;
+                                break;
+                            }
+                        }
+                    }
+
+                    if (is_adj)
+                    {
+                        // basematrix and its data pointer
+                        if (pv.get_row_size(iat1) <= 0 || pv.get_col_size(iat2) <= 0)
+                        {
+                            continue;
+                        }
+                        std::vector<double*> tmp_matrix_ptr;
+                        for (int is = 0; is < GlobalV::NSPIN; ++is)
+                        {
+                            auto* tmp_base_matrix = DM->get_DMR_pointer(is+1)->find_matrix(iat1, iat2, rx2, ry2, rz2);
+                            tmp_matrix_ptr.push_back(tmp_base_matrix->get_pointer());
+                        }
+                        //hamilt::BaseMatrix<double>* tmp_matrix = DM->get_DMR_pointer(1)->find_matrix(iat1, iat2, rx2, ry2, rz2);
+                        //double* tmp_matrix_ptr = tmp_matrix->get_pointer();
+                        for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0)
+                        {
+                            const int T0 = adjs.ntype[ad0];
+                            const int I0 = adjs.natom[ad0];
+                            const int iat = GlobalC::ucell.itia2iat(T0, I0);
+                            if (!iat_recorded && isforce)
+                                adj_iat[ad0] = iat;
+
+                            // mohan add 2010-12-19
+                            if (GlobalC::ucell.infoNL.nproj[T0] == 0)
+                                continue;
+
+                            // const int I0 = GlobalC::GridD.getNatom(ad0);
+                            // const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0);
+                            tau0 = adjs.adjacent_tau[ad0];
+
+                            dtau1 = tau0 - tau1;
+                            dtau2 = tau0 - tau2;
+                            const double distance1 = dtau1.norm2() * pow(GlobalC::ucell.lat0, 2);
+                            const double distance2 = dtau2.norm2() * pow(GlobalC::ucell.lat0, 2);
+
+                            // seems a bug here!! mohan 2011-06-17
+                            rcut1 = pow(GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(),
+                                        2);
+                            rcut2 = pow(GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(),
+                                        2);
+
+                            double r0[3];
+                            double r1[3];
+                            r1[0] = (tau1.x - tau0.x);
+                            r1[1] = (tau1.y - tau0.y);
+                            r1[2] = (tau1.z - tau0.z);
+                            r0[0] = (tau2.x - tau0.x);
+                            r0[1] = (tau2.y - tau0.y);
+                            r0[2] = (tau2.z - tau0.z);
+
+                            if (distance1 >= rcut1 || distance2 >= rcut2)
+                            {
+                                continue;
+                            }
+
+                            const int rx0 = adjs.box[ad0].x;
+                            const int ry0 = adjs.box[ad0].y;
+                            const int rz0 = adjs.box[ad0].z;
+                            key_tuple key1(iat1, -rx0, -ry0, -rz0);
+                            key_tuple key2(iat2, rx2 - rx0, ry2 - ry0, rz2 - rz0);
+
+                            int nnr_inner = 0;
+                            for (int j = 0; j < atom1->nw * GlobalV::NPOL; j++)
+                            {
+                                const int j0 = j / GlobalV::NPOL; // added by zhengdy-soc
+                                const int iw1_all = start1 + j;
+                                const int mu = pv.global2local_row(iw1_all);
+								if (mu < 0)
+								{
+									continue;
+								}
+
+                                for (int k = 0; k < atom2->nw * GlobalV::NPOL; k++)
+                                {
+                                    const int k0 = k / GlobalV::NPOL;
+                                    const int iw2_all = start2 + k;
+                                    const int nu = pv.global2local_col(iw2_all);
+									if (nu < 0)
+									{
+										continue;
+									}
+
+                                    // const Atom* atom0 = &GlobalC::ucell.atoms[T0];
+                                    double nlm[3] = {0, 0, 0};
+                                    std::vector<double> nlm_1 = nlm_tot[iat][key2][iw2_all][0];
+                                    std::vector<std::vector<double>> nlm_2;
+                                    nlm_2.resize(3);
+                                    for (int i = 0; i < 3; i++)
+                                    {
+                                        nlm_2[i] = nlm_tot[iat][key1][iw1_all][i + 1];
+                                    }
+
+                                    assert(nlm_1.size() == nlm_2[0].size());
+
+                                    const int nproj = GlobalC::ucell.infoNL.nproj[T0];
+                                    int ib = 0;
+                                    for (int nb = 0; nb < nproj; nb++)
+                                    {
+                                        const int L0 = GlobalC::ucell.infoNL.Beta[T0].Proj[nb].getL();
+                                        for (int m = 0; m < 2 * L0 + 1; m++)
+                                        {
+                                            for (int ir = 0; ir < 3; ir++)
+                                            {
+                                                nlm[ir] += nlm_2[ir][ib] * nlm_1[ib]
+                                                           * GlobalC::ucell.atoms[T0].ncpp.dion(nb, nb);
+                                            }
+                                            ib += 1;
+                                        }
+                                    }
+                                    assert(ib == nlm_1.size());
+
+                                    double nlm1[3] = {0, 0, 0};
+                                    if (isstress)
+                                    {
+                                        std::vector<double> nlm_1 = nlm_tot[iat][key1][iw1_all][0];
+                                        std::vector<std::vector<double>> nlm_2;
+                                        nlm_2.resize(3);
+                                        for (int i = 0; i < 3; i++)
+                                        {
+                                            nlm_2[i] = nlm_tot[iat][key2][iw2_all][i + 1];
+                                        }
+
+                                        assert(nlm_1.size() == nlm_2[0].size());
+
+                                        const int nproj = GlobalC::ucell.infoNL.nproj[T0];
+                                        int ib = 0;
+                                        for (int nb = 0; nb < nproj; nb++)
+                                        {
+                                            const int L0 = GlobalC::ucell.infoNL.Beta[T0].Proj[nb].getL();
+                                            for (int m = 0; m < 2 * L0 + 1; m++)
+                                            {
+                                                for (int ir = 0; ir < 3; ir++)
+                                                {
+                                                    nlm1[ir] += nlm_2[ir][ib] * nlm_1[ib]
+                                                                * GlobalC::ucell.atoms[T0].ncpp.dion(nb, nb);
+                                                }
+                                                ib += 1;
+                                            }
+                                        }
+                                        assert(ib == nlm_1.size());
+                                    }
+
+                                    /// only one projector for each atom force, but another projector for stress
+                                    force_updated = true;
+                                    // get DMR
+                                    double dm2d1 = 0.0;
+                                    for (int is = 0; is < GlobalV::NSPIN; ++is)
+                                    {  
+                                        dm2d1 += tmp_matrix_ptr[is][nnr_inner];
+                                    }
+                                    double dm2d2 = 2.0 * dm2d1;
+                                    //
+                                    for (int jpol = 0; jpol < 3; jpol++)
+                                    {
+                                        if (isforce)
+                                        {
+                                            adj_fvnl_dbeta[ad0 * 3 + jpol] -= dm2d2 * nlm[jpol];
+                                        }
+                                        if (isstress)
+                                        {
+                                            for (int ipol = jpol; ipol < 3; ipol++)
+                                            {
+                                                local_svnl_dbeta(jpol, ipol)
+                                                    += dm2d1
+                                                        * (nlm[jpol] * r1[ipol] + nlm1[jpol] * r0[ipol]);
+                                            }
+                                        }
+                                    }
+                                    //}
+                                    nnr_inner++;
+                                } // k
+                            }     // j
+                        }         // ad0
+
+                        // outer circle : accumulate nnr
+                        for (int j = 0; j < atom1->nw * GlobalV::NPOL; j++)
+                        {
+                            const int j0 = j / GlobalV::NPOL; // added by zhengdy-soc
+                            const int iw1_all = start1 + j;
+                            const int mu = pv.global2local_row(iw1_all);
+							if (mu < 0)
+							{
+								continue;
+							}
+
+                            // fix a serious bug: atom2[T2] -> atom2
+                            // mohan 2010-12-20
+                            for (int k = 0; k < atom2->nw * GlobalV::NPOL; k++)
+                            {
+                                const int k0 = k / GlobalV::NPOL;
+                                const int iw2_all = start2 + k;
+                                const int nu = pv.global2local_col(iw2_all);
+								if (nu < 0)
+								{
+									continue;
+								}
+								total_nnr++;
+                                nnr++;
+                            }
+                        }
+                        iat_recorded = true;
+                    } // is_adj
+                }     // ad2
+
+                // sum the diff to fvnl_dbeta
+                if (force_updated && isforce)
+                {
+#ifdef _OPENMP
+                    if (num_threads > 1)
+                    {
+                        for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0)
+                        {
+#pragma omp atomic
+                            fvnl_dbeta(adj_iat[ad0], 0) += adj_fvnl_dbeta[ad0 * 3 + 0];
+#pragma omp atomic
+                            fvnl_dbeta(adj_iat[ad0], 1) += adj_fvnl_dbeta[ad0 * 3 + 1];
+#pragma omp atomic
+                            fvnl_dbeta(adj_iat[ad0], 2) += adj_fvnl_dbeta[ad0 * 3 + 2];
+                        }
+                    }
+                    else
+#endif
+                    {
+                        for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0)
+                        {
+                            fvnl_dbeta(adj_iat[ad0], 0) += adj_fvnl_dbeta[ad0 * 3 + 0];
+                            fvnl_dbeta(adj_iat[ad0], 1) += adj_fvnl_dbeta[ad0 * 3 + 1];
+                            fvnl_dbeta(adj_iat[ad0], 2) += adj_fvnl_dbeta[ad0 * 3 + 2];
+                        }
+                    }
+                }
+            } // I1
+        }     // T1
+     
+#ifdef _OPENMP
+        if (isstress)
+        {
+#pragma omp critical(cal_fvnl_dbeta_k_new_reduce)
+            {
+                for (int l = 0; l < 3; l++)
+                {
+                    for (int m = 0; m < 3; m++)
+                    {
+                        svnl_dbeta(l, m) += local_svnl_dbeta(l, m);
+                    }
+                }
+            }
+        }
+    }
+#endif
+
+    assert(total_nnr == pv.nnr);
+
+    if (isstress)
+    {
+        StressTools::stress_fill(GlobalC::ucell.lat0, GlobalC::ucell.omega, svnl_dbeta);
+    }
+
+    ModuleBase::timer::tick("Force_LCAO_k", "cal_fvnl_dbeta_k_new");
+    return;
+}
+
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp
index c648e25fa32..a5fa9a54e3e 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp
@@ -63,7 +63,6 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only)
 	ModuleBase::GlobalFunc::ZEROS(na_each, na_proc);
 	int iat = 0;
 	
-//	std::cout << " in for_2d" << std::endl;
 
 	for (int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1)
 	{
diff --git a/source/module_io/write_HS_R.cpp b/source/module_io/write_HS_R.cpp
index 2fbb56bbca3..81346621f78 100644
--- a/source/module_io/write_HS_R.cpp
+++ b/source/module_io/write_HS_R.cpp
@@ -23,13 +23,13 @@ void ModuleIO::output_HS_R(const int& istep,
     if(GlobalV::NSPIN==1||GlobalV::NSPIN==4)
     {
         // jingan add 2021-6-4, modify 2021-12-2
-        UHM.calculate_HSR_sparse(0, sparse_threshold, kv.nmp, p_ham);
+        UHM.cal_HSR_sparse(0, sparse_threshold, kv.nmp, p_ham);
     }
     else if(GlobalV::NSPIN==2)
     {
         // save HR of current_spin first
-        UHM.calculate_HSR_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, kv.nmp, p_ham);
-        // calculate HR of the other spin
+        UHM.cal_HSR_sparse(GlobalV::CURRENT_SPIN, sparse_threshold, kv.nmp, p_ham);
+        // cal HR of the other spin
         if(GlobalV::VL_IN_H)
         {
             int ik = 0;
@@ -46,7 +46,7 @@ void ModuleIO::output_HS_R(const int& istep,
             p_ham->refresh();
             p_ham->updateHk(ik);
         }
-        UHM.calculate_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);
@@ -70,7 +70,7 @@ void ModuleIO::output_dH_R(const int& istep,
     UHM.GK.allocate_pvdpR();
     if(GlobalV::NSPIN==1||GlobalV::NSPIN==4)
     {
-        UHM.calculate_dH_sparse(0, sparse_threshold);
+        UHM.cal_dH_sparse(0, sparse_threshold);
     }
     else if(GlobalV::NSPIN==2)
     {
@@ -95,7 +95,7 @@ void ModuleIO::output_dH_R(const int& istep,
                     }
                 }
 
-                UHM.calculate_dH_sparse(GlobalV::CURRENT_SPIN, sparse_threshold);
+                UHM.cal_dH_sparse(GlobalV::CURRENT_SPIN, sparse_threshold);
             }
         }
     }
@@ -119,7 +119,7 @@ void ModuleIO::output_S_R(
     ModuleBase::TITLE("ModuleIO","output_S_R");
     ModuleBase::timer::tick("ModuleIO","output_S_R"); 
 
-    UHM.calculate_SR_sparse(sparse_threshold, p_ham);
+    UHM.cal_SR_sparse(sparse_threshold, p_ham);
     ModuleIO::save_SR_sparse(*UHM.LM, sparse_threshold, binary, SR_filename);
     UHM.destroy_all_HSR_sparse();
 
@@ -148,10 +148,10 @@ void ModuleIO::output_T_R(
         sst << GlobalV::global_out_dir << TR_filename;
     }
 
-    UHM.calculate_TR_sparse(sparse_threshold);
+    UHM.cal_TR_sparse(sparse_threshold);
     ModuleIO::save_TR_sparse(istep, *UHM.LM, sparse_threshold, binary, sst.str().c_str());
     UHM.destroy_TR_sparse();
 
     ModuleBase::timer::tick("ModuleIO","output_T_R");
     return;
-}
\ No newline at end of file
+}