diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index ff5fdc1f23..1aea12100a 100644
--- a/source/Makefile.Objects
+++ b/source/Makefile.Objects
@@ -643,6 +643,7 @@ OBJS_LCAO=evolve_elec.o\
       LCAO_init_basis.o\
       setup_exx.o\
       setup_deepks.o\
+      setup_dm.o\
       rho_tau_lcao.o\
       center2_orb.o\
       center2_orb-orb11.o\
diff --git a/source/source_esolver/esolver_dm2rho.cpp b/source/source_esolver/esolver_dm2rho.cpp
index 4c78c5db53..84e6da5141 100644
--- a/source/source_esolver/esolver_dm2rho.cpp
+++ b/source/source_esolver/esolver_dm2rho.cpp
@@ -46,30 +46,22 @@ void ESolver_DM2rho<TK, TR>::runner(UnitCell& ucell, const int istep)
 
     ESolver_KS_LCAO<TK, TR>::before_scf(ucell, istep);
 
-    auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
-
-    if(!estate)
-    {
-        ModuleBase::WARNING_QUIT("ESolver_DM2rho::after_scf","pelec does not exist");
-    }
-
     // file name of DM
     std::string zipname = "output_DM0.npz";
-    elecstate::DensityMatrix<TK, double>* dm = estate->get_DM();
 
     // read DM from file
-    ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1)));
+    ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(1)));
 
     // if nspin=2, need extra reading
     if (PARAM.inp.nspin == 2)
     {
         zipname = "output_DM1.npz";
-        ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2)));
+        ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(2)));
     }
 
     // it's dangerous to design psiToRho function like this, mohan note 20251024
     // this->pelec->psiToRho(*this->psi);
-    LCAO_domain::dm2rho(estate->DM->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
+    LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
 
     int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1;
 
diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp
index 814c675b8c..082271ee2e 100644
--- a/source/source_esolver/esolver_double_xc.cpp
+++ b/source/source_esolver/esolver_double_xc.cpp
@@ -84,7 +84,7 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
     }
 
     // 6) initialize the density matrix
-    dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->init_DM(&this->kv, &(this->pv), PARAM.inp.nspin);
+    this->dmat_base.allocate_dm(&this->kv, &this->pv, PARAM.inp.nspin);
 
     // 10) inititlize the charge density
     this->chr_base.allocate(PARAM.inp.nspin);
@@ -138,8 +138,6 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
     if (this->p_hamilt_base == nullptr)
     {
-        elecstate::DensityMatrix<TK, double>* DM = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->get_DM();
-
         this->p_hamilt_base = new hamilt::HamiltLCAO<TK, TR>(
             ucell,
             this->gd,
@@ -148,7 +146,7 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
             this->kv,
             this->two_center_bundle_,
             this->orb_,
-            DM,
+            this->dmat_base.dm,
 			this->deepks,
 			istep,
 			this->exx_nao);
@@ -159,13 +157,11 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); 
 
     // DMR should be same size with Hamiltonian(R)
-    dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)
-    ->get_DM()
-    ->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt_base)->getHR()));
+    this->dmat_base.dm->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt_base)->getHR()));
 
     if (istep > 0)
     {
-        dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->get_DM()->cal_DMR();
+        this->dmat_base.dm->cal_DMR();
     }
 
     ModuleBase::timer::tick("ESolver_DoubleXC", "before_scf");
@@ -226,23 +222,23 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
         std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&this->deepks.ld, [](LCAO_Deepks<TK>*) {});
         LCAO_Deepks_Interface<TK, TR> deepks_interface(ld_shared_ptr);
 
-        deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
-                                            this->kv.get_nks(),
-                                            ucell.nat,
-                                            PARAM.globalv.nlocal,
-                                            this->pelec->ekb,
-                                            this->kv.kvec_d,
-                                            ucell,
-                                            this->orb_,
-                                            this->gd,
-                                            &(this->pv),
-                                            *(this->psi),
-                                            dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
-                                            p_ham_deepks,
-                                            iter,
-                                            conv_esolver,
-                                            GlobalV::MY_RANK,
-                                            GlobalV::ofs_running);
+		deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
+				this->kv.get_nks(),
+				ucell.nat,
+				PARAM.globalv.nlocal,
+				this->pelec->ekb,
+				this->kv.kvec_d,
+				ucell,
+				this->orb_,
+				this->gd,
+				&(this->pv),
+				*(this->psi),
+				this->dmat.dm,
+				p_ham_deepks,
+				iter,
+				conv_esolver,
+				GlobalV::MY_RANK,
+				GlobalV::ofs_running);
 #endif
                                             
         // restore to density after charge mixing
@@ -352,9 +348,12 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
             auto _pes_lcao = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
             for (int ik = 0; ik < nks; ik++)
             {
-                _pes_lcao_base->get_DM()->set_DMK_pointer(ik, _pes_lcao->get_DM()->get_DMK_pointer(ik));
+// mohan update 2025-11-03
+                this->dmat_base.dm->set_DMK_pointer(ik, this->dmat.dm->get_DMK_pointer(ik));
+//                _pes_lcao_base->get_DM()->set_DMK_pointer(ik, _pes_lcao->get_DM()->get_DMK_pointer(ik));
             }
-            _pes_lcao_base->get_DM()->cal_DMR();
+            this->dmat_base.dm->cal_DMR();
+//            _pes_lcao_base->get_DM()->cal_DMR();
             _pes_lcao_base->ekb = _pes_lcao->ekb;
             _pes_lcao_base->wg = _pes_lcao->wg;          
         }        
@@ -386,6 +385,7 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                        this->gd,
                        this->pv,
                        this->pelec_base,
+                       this->dmat_base, // mohan add 2025-11-03
                        this->psi,
                        this->two_center_bundle_,
                        this->orb_,
diff --git a/source/source_esolver/esolver_double_xc.h b/source/source_esolver/esolver_double_xc.h
index 591411c82d..bff9e28bc7 100644
--- a/source/source_esolver/esolver_double_xc.h
+++ b/source/source_esolver/esolver_double_xc.h
@@ -32,6 +32,9 @@ class ESolver_DoubleXC : public ESolver_KS_LCAO<TK, TR>
     //! Electronic states
     elecstate::ElecState* pelec_base = nullptr;
 
+    //! Density Matrix, mohan add 2025-11-03
+    LCAO_domain::Setup_DM<TK> dmat_base;
+
     //! Electorn charge density
     Charge chr_base;
 };
diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp
index ef2258de36..cca74c26a5 100644
--- a/source/source_esolver/esolver_fp.cpp
+++ b/source/source_esolver/esolver_fp.cpp
@@ -54,9 +54,7 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp)
 
     // write geometry file
     ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif",
-                               ucell,
-                               "# Generated by ABACUS ModuleIO::CifParser",
-                               "data_?");
+      ucell, "# Generated by ABACUS ModuleIO::CifParser", "data_?");
 
     // init charge extrapolation
     this->CE.Init_CE(inp.nspin, ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap);
diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp
index f85fd53f43..2818152b22 100644
--- a/source/source_esolver/esolver_ks.cpp
+++ b/source/source_esolver/esolver_ks.cpp
@@ -45,7 +45,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 {
     ModuleBase::TITLE("ESolver_KS", "before_all_runners");
 
-    //! 1) init "before_all_runniers" in ESolver_FP
+    //! 1) setup "before_all_runniers" in ESolver_FP
     ESolver_FP::before_all_runners(ucell, inp);
     
     //! 2) setup some parameters
@@ -67,7 +67,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
 
-    //! 4) setup Exc for the first element '0', because all elements have same exc 
+    //! 4) setup Exc for the first element '0' (all elements have same exc) 
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
     
     //! 5) setup the charge mixing parameters
@@ -84,7 +84,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
     }
 
-    //! 7) Setup the k points according to symmetry.
+    //! 7) setup k points in the Brillouin zone according to symmetry.
     this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
@@ -99,19 +99,12 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 			this->pw_rhod->nplane, this->pw_rhod->nrxx, pw_big->nbz, pw_big->bz);
 
     //! 11) calculate the structure factor
-    this->sf.setup_structure_factor(&ucell, Pgrid, this->pw_rhod);
+    this->sf.setup(&ucell, Pgrid, this->pw_rhod);
 }
 
 template <typename T, typename Device>
 void ESolver_KS<T, Device>::hamilt2rho_single(UnitCell& ucell, const int istep, const int iter, const double ethr)
-{
-    ModuleBase::timer::tick(this->classname, "hamilt2rho_single");
-    // Temporarily, before HSolver is constructed, it should be overrided by
-    // LCAO, PW, SDFT and TDDFT.
-    // After HSolver is constructed, LCAO, PW, SDFT should delete their own
-    // hamilt2rho_single() and use:
-    ModuleBase::timer::tick(this->classname, "hamilt2rho_single");
-}
+{}
 
 template <typename T, typename Device>
 void ESolver_KS<T, Device>::hamilt2rho(UnitCell& ucell, const int istep, const int iter, const double ethr)
diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp
index 47382550b4..b3f83c42b8 100644
--- a/source/source_esolver/esolver_ks_lcao.cpp
+++ b/source/source_esolver/esolver_ks_lcao.cpp
@@ -91,7 +91,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
 
 
     // 7) init DMK, but DMR is constructed in before_scf()
-    dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->init_DM(&this->kv, &(this->pv), inp.nspin);
+    this->dmat.allocate_dm(&this->kv, &this->pv, inp.nspin);
 
     // 8) init exact exchange calculations
     this->exx_nao.before_runner(ucell, this->kv, this->orb_, this->pv, inp);
@@ -153,12 +153,6 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     //! 1) call before_scf() of ESolver_KS.
     ESolver_KS<TK>::before_scf(ucell, istep);
 
-    auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
-    if(!estate)
-    {
-        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf","pelec does not exist");
-    }
-
     //! 2) find search radius
     double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
       PARAM.inp.out_level, orb_.get_rcutmax_Phi(), ucell.infoNL.get_rcutmax_Beta(),
@@ -193,11 +187,9 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
     if (this->p_hamilt == nullptr)
     {
-        elecstate::DensityMatrix<TK, double>* DM = estate->get_DM();
-
         this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
             ucell, this->gd, &this->pv, this->pelec->pot, this->kv,
-            two_center_bundle_, orb_, DM, this->deepks, istep, exx_nao);
+            two_center_bundle_, orb_, this->dmat.dm, this->deepks, istep, exx_nao);
     }
 
     // 9) for each ionic step, the overlap <phi|alpha> must be rebuilt
@@ -210,7 +202,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
         sc.init_sc(PARAM.inp.sc_thr, PARAM.inp.nsc, PARAM.inp.nsc_min, PARAM.inp.alpha_trial,
                    PARAM.inp.sccut, PARAM.inp.sc_drop_thr, ucell, &(this->pv),
-                   PARAM.inp.nspin, this->kv, this->p_hamilt, this->psi, this->pelec);
+                   PARAM.inp.nspin, this->kv, this->p_hamilt, this->psi, this->dmat.dm, this->pelec);
     }
 
     // 11) set xc type before the first cal of xc in pelec->init_scf, Peize Lin add 2016-12-03
@@ -219,17 +211,16 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     // 12) init_scf, should be before_scf? mohan add 2025-03-10
     this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
 
-    // 13) initalize DMR
-    // DMR should be same size with Hamiltonian(R)
+    // 13) initalize DM(R), which has the same size with Hamiltonian(R)
     auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
     if(!hamilt_lcao)
     {
         ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf","p_hamilt does not exist");
     }
-    estate->get_DM()->init_DMR(*hamilt_lcao->getHR());
+    this->dmat.dm->init_DMR(*hamilt_lcao->getHR());
 
 #ifdef __MLALGO
-    // 14) initialize DMR of DeePKS
+    // 14) initialize DM2(R) of DeePKS, the DM2(R) is different from DM(R)
     this->deepks.ld.init_DMR(ucell, orb_, this->pv, this->gd);
 #endif
 
@@ -238,7 +229,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
     if (istep > 0)
     {
-        estate->get_DM()->cal_DMR();
+        this->dmat.dm->cal_DMR();
     }
 
     // 16) the electron charge density should be symmetrized,
@@ -265,8 +256,6 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
 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;
 }
 
@@ -282,7 +271,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
 
     fsl.getForceStress(ucell, PARAM.inp.cal_force, PARAM.inp.cal_stress, 
                        PARAM.inp.test_force, PARAM.inp.test_stress,
-                       this->gd, this->pv, this->pelec, this->psi,
+                       this->gd, this->pv, this->pelec, this->dmat, this->psi,
                        two_center_bundle_, orb_, force, this->scs,
                        this->locpp, this->sf, this->kv,
                        this->pw_rho, this->solvent, this->deepks,
@@ -302,15 +291,13 @@ void ESolver_KS_LCAO<TK, TR>::cal_stress(UnitCell& ucell, ModuleBase::matrix& st
     ModuleBase::TITLE("ESolver_KS_LCAO", "cal_stress");
     ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress");
 
-    // if the users do not want to calculate forces but want stress,
-    // we call cal_force
     if (!this->have_force)
     {
         ModuleBase::matrix fcs;
         this->cal_force(ucell, fcs);
     }
 
-    // the 'scs' stress has already been calculated in 'cal_force'
+    // the stress has been calculated in 'cal_force'
     stress = this->scs;
     this->have_force = false;
 
@@ -327,21 +314,14 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
 
     const int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1;
 
-	auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
     auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
-
-	if(!estate)
-	{
-		ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_all_runners","pelec does not exist");
-	}
-
 	if(!hamilt_lcao)
 	{
 		ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_all_runners","p_hamilt does not exist");
 	}
 
 	ModuleIO::ctrl_runner_lcao<TK, TR>(ucell,
-		  PARAM.inp, this->kv, estate, this->pv, this->Pgrid, 
+		  PARAM.inp, this->kv, this->pelec, this->dmat, this->pv, this->Pgrid, 
 		  this->gd, this->psi, this->chr, hamilt_lcao,
           this->two_center_bundle_,
           this->orb_, this->pw_rho, this->pw_rhod,
@@ -358,19 +338,10 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     // call iter_init() of ESolver_KS
     ESolver_KS<TK>::iter_init(ucell, istep, iter);
 
-    // cast pointers
-	auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
-	if(!estate)
-	{
-		ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::iter_init","pelec does not exist");
-	}
-
-	elecstate::DensityMatrix<TK, double>* dm = estate->get_DM();
-
-    module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp); 
+    module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, this->dmat.dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp); 
 
     // mohan update 2012-06-05
-    estate->f_en.deband_harris = estate->cal_delta_eband(ucell);
+    this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell);
 
     if (istep == 0 && PARAM.inp.init_wfc == "file")
 	{
@@ -384,7 +355,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
                   this->exx_nao.exd->two_level_step : this->exx_nao.exc->two_level_step;
 		}
 #endif
-		elecstate::init_dm<TK>(ucell, estate, this->psi, this->chr, iter, exx_two_level_step);
+		elecstate::init_dm<TK>(ucell, this->pelec, this->dmat, this->psi, this->chr, iter, exx_two_level_step);
 	}
 
 #ifdef __EXX
@@ -393,11 +364,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     {
         if (GlobalC::exx_info.info_ri.real_number)
         {
-            this->exx_nao.exd->exx_eachiterinit(istep, ucell, *dm, this->kv, iter);
+            this->exx_nao.exd->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, iter);
         }
         else
         {
-            this->exx_nao.exc->exx_eachiterinit(istep, ucell, *dm, this->kv, iter);
+            this->exx_nao.exc->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, iter);
         }
     }
 #endif
@@ -406,7 +377,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     {
         if (istep != 0 || iter != 1)
         {
-            GlobalC::dftu.set_dmr(dm);
+            GlobalC::dftu.set_dmr(this->dmat.dm);
         }
         // Calculate U and J if Yukawa potential is used
         GlobalC::dftu.cal_slater_UJ(ucell, this->chr.rho, this->pw_rho->nrxx);
@@ -432,7 +403,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     // save density matrix DMR for mixing
     if (PARAM.inp.mixing_restart > 0 && PARAM.inp.mixing_dmr && this->p_chgmix->mixing_restart_count > 0)
     {
-        dm->save_DMR();
+        this->dmat.dm->save_DMR();
     }
 }
 
@@ -470,7 +441,8 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
     if (!skip_solve)
     {
         hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
-        hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, this->chr, PARAM.inp.nspin, skip_charge);
+        hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, *this->dmat.dm, 
+          this->chr, PARAM.inp.nspin, skip_charge);
     }
 
     // 4) EXX
@@ -505,20 +477,14 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
 {
     ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish");
 
-    auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
     auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
 
-    if(!estate)
-    {
-        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::iter_finish","pelec does not exist");
-    }
-
     if(!hamilt_lcao)
     {
         ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::iter_finish","p_hamilt does not exist");
     }
 
-	const std::vector<std::vector<TK>>& dm_vec = estate->get_DM()->get_DMK_vector();
+	const std::vector<std::vector<TK>>& dm_vec = this->dmat.dm->get_DMK_vector();
 
     // 1) calculate the local occupation number matrix and energy correction in DFT+U
     if (PARAM.inp.dft_plus_u)
@@ -559,8 +525,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     {
         if (PARAM.inp.mixing_restart > 0 && this->p_chgmix->mixing_restart_count > 0 && PARAM.inp.mixing_dmr)
         {
-            elecstate::DensityMatrix<TK, double>* dm = estate->get_DM();
-            this->p_chgmix->mix_dmr(dm);
+            this->p_chgmix->mix_dmr(this->dmat.dm);
         }
     }
 
@@ -571,11 +536,10 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     }
 
     // control the output related to the finished iteration
-    ModuleIO::ctrl_iter_lcao<TK, TR>(ucell, PARAM.inp, this->kv, estate,
+    ModuleIO::ctrl_iter_lcao<TK, TR>(ucell, PARAM.inp, this->kv, this->pelec, *this->dmat.dm,
       this->pv, this->gd, this->psi, this->chr, this->p_chgmix, 
       hamilt_lcao, this->orb_, this->deepks, 
       this->exx_nao, iter, istep, conv_esolver, this->scf_ene_thr);
-
 }
 
 template <typename TK, typename TR>
@@ -584,14 +548,8 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
     ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf");
     ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf");
 
-    auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
     auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
 
-    if(!estate)
-    {
-        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_scf","pelec does not exist");
-    }
-
     if(!hamilt_lcao)
     {
         ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_scf","p_hamilt does not exist");
@@ -599,24 +557,19 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
 
     if (PARAM.inp.out_elf[0] > 0)
 	{
-		LCAO_domain::dm2tau(estate->DM->get_DMR_vector(), PARAM.inp.nspin, estate->charge);
+		LCAO_domain::dm2tau(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, this->pelec->charge);
 	}
 
     //! 1) call after_scf() of ESolver_KS
     ESolver_KS<TK>::after_scf(ucell, istep, conv_esolver);
 
-
     //! 2) output of lcao every few ionic steps
     ModuleIO::ctrl_scf_lcao<TK, TR>(ucell,
-            PARAM.inp, this->kv, estate, this->pv,
-            this->gd, this->psi, hamilt_lcao,
-            this->two_center_bundle_,
-            this->orb_, this->pw_wfc, this->pw_rho,
-            this->pw_big, this->sf,
+            PARAM.inp, this->kv, this->pelec, this->dmat.dm, this->pv,
+            this->gd, this->psi, hamilt_lcao, this->two_center_bundle_,
+            this->orb_, this->pw_wfc, this->pw_rho, this->pw_big, this->sf,
             this->rdmft_solver, this->deepks, this->exx_nao, 
-            this->conv_esolver, this->scf_nmax_flag,
-            istep);
-
+            this->conv_esolver, this->scf_nmax_flag, istep);
 
     //! 3) Clean up RA, which is used to serach for adjacent atoms
     if (!PARAM.inp.cal_force && !PARAM.inp.cal_stress)
diff --git a/source/source_esolver/esolver_ks_lcao.h b/source/source_esolver/esolver_ks_lcao.h
index 16fb6aa8df..2bcf0c672c 100644
--- a/source/source_esolver/esolver_ks_lcao.h
+++ b/source/source_esolver/esolver_ks_lcao.h
@@ -9,6 +9,7 @@
 #include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251008
 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008
 #include "source_lcao/module_rdmft/rdmft.h" // rdmft
+#include "source_lcao/setup_dm.h" // mohan add 2025-10-30
 
 #include <memory>
 
@@ -71,6 +72,9 @@ class ESolver_KS_LCAO : public ESolver_KS<TK>
     //! NAO orbitals: two-center integrations
     TwoCenterBundle two_center_bundle_;
 
+    //! Add density matrix class, mohan add 2025-10-30
+    LCAO_domain::Setup_DM<TK> dmat;
+
     //! For RDMFT calculations, added by jghan, 2024-03-16 
     rdmft::RDMFT<TK, TR> rdmft_solver;
 
diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp
index 6e864b4e4a..d608082ddd 100644
--- a/source/source_esolver/esolver_ks_lcao_tddft.cpp
+++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp
@@ -264,7 +264,8 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
         {
             bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
             hsolver::HSolverLCAO<std::complex<double>> hsolver_lcao_obj(&this->pv, PARAM.inp.ks_solver);
-            hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, this->chr, PARAM.inp.nspin, skip_charge);
+			hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, *this->dmat.dm, 
+					this->chr, PARAM.inp.nspin, skip_charge);
         }
     }
 
@@ -317,7 +318,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::iter_finish(UnitCell& ucell,
     if (conv_esolver && estep == estep_max - 1 && istep >= (PARAM.inp.init_wfc == "file" ? 0 : 1)
         && PARAM.inp.td_edm == 0)
     {
-        elecstate::cal_edm_tddft(this->pv, this->pelec, this->kv, this->p_hamilt);
+        elecstate::cal_edm_tddft(this->pv, this->dmat, this->kv, this->p_hamilt);
     }
 }
 
@@ -431,9 +432,9 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::store_h_s_psi(UnitCell& ucell,
                                     1,
                                     this->Sk_laststep.template data<std::complex<double>>() + ik * len_HS_ik,
                                     1);
-            }
-        }
-    }
+            } // end use_tensor
+        } // end ik
+    }// conv_esolver
 }
 
 template <typename TR, typename Device>
@@ -481,20 +482,18 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::weight_dm_rho(const UnitCell& ucell)
     // Calculate Eband energy
     elecstate::calEBand(this->pelec->ekb, this->pelec->wg, this->pelec->f_en);
 
-    // Calculate the density matrix
-    auto _pes = dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec);
-    elecstate::cal_dm_psi(_pes->DM->get_paraV_pointer(), _pes->wg, this->psi[0], *(_pes->DM));
-    if (PARAM.inp.td_stype == 2)
+    elecstate::cal_dm_psi(this->dmat.dm->get_paraV_pointer(), this->pelec->wg, this->psi[0], *this->dmat.dm);
+    if(PARAM.inp.td_stype == 2)
     {
-        _pes->DM->cal_DMR_td(ucell, TD_info::cart_At);
+        this->dmat.dm->cal_DMR_td(ucell, TD_info::cart_At);
     }
     else
     {
-        _pes->DM->cal_DMR();
+        this->dmat.dm->cal_DMR();
     }
 
-    // Get the real-space charge density, mohan add 2025-10-24
-    LCAO_domain::dm2rho(_pes->DM->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
+    // get the real-space charge density, mohan add 2025-10-24
+    LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
 }
 
 template class ESolver_KS_LCAO_TDDFT<double, base_device::DEVICE_CPU>;
diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp
index 6785c64a1f..e3fad24322 100644
--- a/source/source_esolver/esolver_of.cpp
+++ b/source/source_esolver/esolver_of.cpp
@@ -105,7 +105,7 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp)
                         pw_big->nbz,
                         pw_big->bz); // mohan add 2010-07-22, update 2011-05-04
     // Calculate Structure factor
-    sf.setup_structure_factor(&ucell, Pgrid, pw_rho);
+    sf.setup(&ucell, Pgrid, pw_rho);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS");
 
     // initialize local pseudopotential
diff --git a/source/source_esolver/lcao_others.cpp b/source/source_esolver/lcao_others.cpp
index a31dc28a9c..9c413836d4 100644
--- a/source/source_esolver/lcao_others.cpp
+++ b/source/source_esolver/lcao_others.cpp
@@ -173,7 +173,6 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     }
     if (this->p_hamilt == nullptr)
     {
-        elecstate::DensityMatrix<TK, double>* DM = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
         this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
             ucell,
             this->gd,
@@ -182,7 +181,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
             this->kv,
             two_center_bundle_,
             orb_,
-            DM,
+            this->dmat.dm,
 			this->deepks,
 			istep,
 			this->exx_nao);
@@ -208,6 +207,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
                    this->kv,
                    this->p_hamilt,
                    this->psi,
+                   this->dmat.dm,
                    this->pelec);
     }
 
diff --git a/source/source_estate/elecstate_lcao.cpp b/source/source_estate/elecstate_lcao.cpp
index b0518d817e..14562fca87 100644
--- a/source/source_estate/elecstate_lcao.cpp
+++ b/source/source_estate/elecstate_lcao.cpp
@@ -15,12 +15,6 @@
 namespace elecstate
 {
 
-template <typename TK>
-void ElecStateLCAO<TK>::init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin)
-{
-    const int nspin_dm = nspin == 2 ? 2 : 1;
-    this->DM = new DensityMatrix<TK, double>(paraV, nspin_dm, kv->kvec_d, kv->get_nks() / nspin_dm);
-}
 
 template <>
 double ElecStateLCAO<double>::get_spin_constrain_energy()
diff --git a/source/source_estate/elecstate_lcao.h b/source/source_estate/elecstate_lcao.h
index e9bb3077a5..5d2d9ac7f6 100644
--- a/source/source_estate/elecstate_lcao.h
+++ b/source/source_estate/elecstate_lcao.h
@@ -27,21 +27,11 @@ class ElecStateLCAO : public ElecState
 
     virtual ~ElecStateLCAO()
     {
-        if (this->DM != nullptr)
-        {
-            delete this->DM;
-        }
     }
 
     // update charge density for next scf step
     // void getNewRho() override;
 
-    // initial density matrix
-    void init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin);
-    DensityMatrix<TK, double>* get_DM() const
-    {
-        return const_cast<DensityMatrix<TK, double>*>(this->DM);
-    }
     static int out_wfc_lcao;
     static bool need_psi_grid;
 
@@ -59,8 +49,6 @@ class ElecStateLCAO : public ElecState
     void dm2rho(std::vector<TK*> pexsi_DM, std::vector<TK*> pexsi_EDM);
 #endif
 
-    DensityMatrix<TK, double>* DM = nullptr;
-
 };
 
 template <typename TK>
diff --git a/source/source_estate/module_charge/charge_extra.cpp b/source/source_estate/module_charge/charge_extra.cpp
index 345a34b723..dd7bea8344 100644
--- a/source/source_estate/module_charge/charge_extra.cpp
+++ b/source/source_estate/module_charge/charge_extra.cpp
@@ -107,7 +107,7 @@ void Charge_Extra::extrapolate_charge(
     rho_extr = std::min(istep, pot_order);
     if(rho_extr == 0)
     {
-        sf->setup_structure_factor(&ucell, *Pgrid, chr->rhopw);
+        sf->setup(&ucell, *Pgrid, chr->rhopw);
         ofs_running << " charge density from previous step !" << std::endl;
         return;
     }
@@ -169,7 +169,7 @@ void Charge_Extra::extrapolate_charge(
         }
     }
 
-    sf->setup_structure_factor(&ucell, *Pgrid, chr->rhopw);
+    sf->setup(&ucell, *Pgrid, chr->rhopw);
     double** rho_atom = new double*[this->nspin];
     for (int is = 0; is < this->nspin; is++)
     {
@@ -314,4 +314,4 @@ void Charge_Extra::update_delta_rho(const UnitCell& ucell, const Charge* chr, co
     }
     delete[] rho_atom;
     return;
-}
\ No newline at end of file
+}
diff --git a/source/source_estate/module_dm/cal_edm_tddft.cpp b/source/source_estate/module_dm/cal_edm_tddft.cpp
index 9a90ee4cf7..809d76e1e4 100644
--- a/source/source_estate/module_dm/cal_edm_tddft.cpp
+++ b/source/source_estate/module_dm/cal_edm_tddft.cpp
@@ -2,11 +2,13 @@
 
 #include "source_base/module_external/lapack_connector.h"
 #include "source_base/module_external/scalapack_connector.h"
+
+#include "source_io/module_parameter/parameter.h" // use PARAM.globalv
 namespace elecstate
 {
 // use the original formula (Hamiltonian matrix) to calculate energy density matrix
 void cal_edm_tddft(Parallel_Orbitals& pv,
-                   elecstate::ElecState* pelec,
+                   LCAO_domain::Setup_DM<std::complex<double>> &dmat,
                    K_Vectors& kv,
                    hamilt::Hamilt<std::complex<double>>* p_hamilt)
 {
@@ -14,15 +16,13 @@ void cal_edm_tddft(Parallel_Orbitals& pv,
     const int nlocal = PARAM.globalv.nlocal;
     assert(nlocal >= 0);
 
-    auto _pelec = dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(pelec);
-
-    _pelec->get_DM()->EDMK.resize(kv.get_nks());
+    dmat.dm->EDMK.resize(kv.get_nks());
 
     for (int ik = 0; ik < kv.get_nks(); ++ik)
     {
         p_hamilt->updateHk(ik);
-        std::complex<double>* tmp_dmk = _pelec->get_DM()->get_DMK_pointer(ik);
-        ModuleBase::ComplexMatrix& tmp_edmk = _pelec->get_DM()->EDMK[ik];
+        std::complex<double>* tmp_dmk = dmat.dm->get_DMK_pointer(ik);
+        ModuleBase::ComplexMatrix& tmp_edmk = dmat.dm->EDMK[ik];
 
 #ifdef __MPI
 
diff --git a/source/source_estate/module_dm/cal_edm_tddft.h b/source/source_estate/module_dm/cal_edm_tddft.h
index 49de0550b5..5ceec1bdaa 100644
--- a/source/source_estate/module_dm/cal_edm_tddft.h
+++ b/source/source_estate/module_dm/cal_edm_tddft.h
@@ -3,13 +3,13 @@
 
 #include "source_basis/module_ao/parallel_orbitals.h"
 #include "source_cell/klist.h"
-#include "source_estate/elecstate_lcao.h"
 #include "source_hamilt/hamilt.h"
+#include "source_lcao/setup_dm.h"
 
 namespace elecstate
 {
 void cal_edm_tddft(Parallel_Orbitals& pv,
-                   elecstate::ElecState* pelec,
+                   LCAO_domain::Setup_DM<std::complex<double>> &dmat,
                    K_Vectors& kv,
                    hamilt::Hamilt<std::complex<double>>* p_hamilt);
 } // namespace elecstate
diff --git a/source/source_estate/module_dm/init_dm.cpp b/source/source_estate/module_dm/init_dm.cpp
index 1913914226..669360ee1e 100644
--- a/source/source_estate/module_dm/init_dm.cpp
+++ b/source/source_estate/module_dm/init_dm.cpp
@@ -5,7 +5,8 @@
 
 template <typename TK>
 void elecstate::init_dm(UnitCell& ucell,
-		elecstate::ElecStateLCAO<TK>* pelec,
+		elecstate::ElecState* pelec,
+        LCAO_domain::Setup_DM<TK> &dmat,
         psi::Psi<TK>* psi,
 		Charge &chr,
         const int iter,
@@ -30,10 +31,10 @@ void elecstate::init_dm(UnitCell& ucell,
 				pelec->skip_weights);
 
 		elecstate::calEBand(pelec->ekb, pelec->wg, pelec->f_en);
-		elecstate::cal_dm_psi(pelec->DM->get_paraV_pointer(), pelec->wg, *psi, *(pelec->DM));
-		pelec->DM->cal_DMR();
+		elecstate::cal_dm_psi(dmat.dm->get_paraV_pointer(), pelec->wg, *psi, *dmat.dm);
+		dmat.dm->cal_DMR();
 
-		pelec->psiToRho(*psi);
+//		pelec->psiToRho(*psi); // I found this sentence is useless, mohan add 2025-11-04
 		pelec->skip_weights = false;
 
 		elecstate::cal_ux(ucell);
@@ -50,14 +51,16 @@ void elecstate::init_dm(UnitCell& ucell,
 
 
 template void elecstate::init_dm<double>(UnitCell& ucell,
-		elecstate::ElecStateLCAO<double>* pelec,
+		elecstate::ElecState* pelec,
+        LCAO_domain::Setup_DM<double> &dmat,
         psi::Psi<double>* psi,
 		Charge &chr,
         const int iter,
         const int exx_two_level_step);
 
 template void elecstate::init_dm<std::complex<double>>(UnitCell& ucell,
-		elecstate::ElecStateLCAO<std::complex<double>>* pelec,
+		elecstate::ElecState* pelec,
+        LCAO_domain::Setup_DM<std::complex<double>> &dmat,
         psi::Psi<std::complex<double>>* psi,
 		Charge &chr,
         const int iter,
diff --git a/source/source_estate/module_dm/init_dm.h b/source/source_estate/module_dm/init_dm.h
index 38b6263eca..2fd969638d 100644
--- a/source/source_estate/module_dm/init_dm.h
+++ b/source/source_estate/module_dm/init_dm.h
@@ -2,16 +2,18 @@
 #define INIT_DM_H
 
 #include "source_cell/unitcell.h" // use unitcell
-#include "source_estate/elecstate_lcao.h"// use ElecStateLCAO
+#include "source_estate/elecstate.h"// use ElecState
 #include "source_psi/psi.h" // use electronic wave functions
 #include "source_estate/module_charge/charge.h" // use charge
+#include "source_lcao/setup_dm.h" // define Setup_DM
 
 namespace elecstate
 {
 
 template <typename TK> 
 void init_dm(UnitCell& ucell,
-		ElecStateLCAO<TK>* pelec,
+		ElecState* pelec,
+        LCAO_domain::Setup_DM<TK> &dmat,
         psi::Psi<TK>* psi,
 		Charge &chr,
         const int iter,
diff --git a/source/source_estate/test/charge_extra_test.cpp b/source/source_estate/test/charge_extra_test.cpp
index 255c2fbe78..c2951c1b5a 100644
--- a/source/source_estate/test/charge_extra_test.cpp
+++ b/source/source_estate/test/charge_extra_test.cpp
@@ -89,7 +89,7 @@ Structure_Factor::Structure_Factor()
 Structure_Factor::~Structure_Factor()
 {
 }
-void Structure_Factor::setup_structure_factor(const UnitCell*, const Parallel_Grid&, const ModulePW::PW_Basis*)
+void Structure_Factor::setup(const UnitCell*, const Parallel_Grid&, const ModulePW::PW_Basis*)
 {
 }
 
diff --git a/source/source_hamilt/module_surchem/cal_pseudo.cpp b/source/source_hamilt/module_surchem/cal_pseudo.cpp
index 1e0a64c5b5..9f707654d0 100644
--- a/source/source_hamilt/module_surchem/cal_pseudo.cpp
+++ b/source/source_hamilt/module_surchem/cal_pseudo.cpp
@@ -9,7 +9,8 @@ void surchem::gauss_charge(const UnitCell& cell,
                            std::complex<double>* N,
                            Structure_Factor* sf)
 {
-    sf->setup_structure_factor(&cell, pgrid, rho_basis); // call strucFac(ntype,ngmc)
+    sf->setup(&cell, pgrid, rho_basis); // this is strange, should be removed to other places, mohan add 2025-11-04
+
     const int ig0 = rho_basis->ig_gge0; // G=0 index
     for (int it = 0; it < cell.ntype; it++)
     {
diff --git a/source/source_hsolver/hsolver_lcao.cpp b/source/source_hsolver/hsolver_lcao.cpp
index bee55de5c5..87eec7475e 100644
--- a/source/source_hsolver/hsolver_lcao.cpp
+++ b/source/source_hsolver/hsolver_lcao.cpp
@@ -39,11 +39,12 @@
 namespace hsolver
 {
 
-template <typename T, typename Device>
-void HSolverLCAO<T, Device>::solve(hamilt::Hamilt<T>* pHamilt,
-                                   psi::Psi<T>& psi,
-                                   elecstate::ElecState* pes,
-                                   Charge &chr,
+template <typename TK, typename Device>
+void HSolverLCAO<TK, Device>::solve(hamilt::Hamilt<TK>* pHamilt,
+                                   psi::Psi<TK>& psi,
+								   elecstate::ElecState* pes,
+								   elecstate::DensityMatrix<TK, double>& dm, // mohan add 2025-11-03
+								   Charge &chr,
                                    const int nspin,
                                    const bool skip_charge)
 {
@@ -94,15 +95,14 @@ void HSolverLCAO<T, Device>::solve(hamilt::Hamilt<T>* pHamilt,
                                      pes->nelec_spin,
                                      pes->skip_weights);
 
-        auto _pes_lcao = dynamic_cast<elecstate::ElecStateLCAO<T>*>(pes);
-        elecstate::calEBand(_pes_lcao->ekb, _pes_lcao->wg, _pes_lcao->f_en);
-        elecstate::cal_dm_psi(_pes_lcao->DM->get_paraV_pointer(), _pes_lcao->wg, psi, *(_pes_lcao->DM));
-        _pes_lcao->DM->cal_DMR();
+        elecstate::calEBand(pes->ekb, pes->wg, pes->f_en);
+        elecstate::cal_dm_psi(dm.get_paraV_pointer(), pes->wg, psi, dm);
+        dm.cal_DMR();
 
         if (!skip_charge)
         {
             // compute charge density from density matrix, mohan update 20251024
-            LCAO_domain::dm2rho(_pes_lcao->DM->get_DMR_vector(), nspin, &chr);
+            LCAO_domain::dm2rho(dm.get_DMR_vector(), nspin, &chr);
         }
         else
         {
@@ -124,7 +124,7 @@ void HSolverLCAO<T, Device>::solve(hamilt::Hamilt<T>* pHamilt,
         auto _pes = dynamic_cast<elecstate::ElecStateLCAO<T>*>(pes);
         pes->f_en.eband = pe.totalFreeEnergy;
         // maybe eferm could be dealt with in the future
-        _pes->dm2rho(pe.DM, pe.EDM);
+        _pes->dm2rho(dm, pe.EDM);
 #endif
     }
 
diff --git a/source/source_hsolver/hsolver_lcao.h b/source/source_hsolver/hsolver_lcao.h
index 67ac1fe32a..eebab34208 100644
--- a/source/source_hsolver/hsolver_lcao.h
+++ b/source/source_hsolver/hsolver_lcao.h
@@ -6,31 +6,33 @@
 #include "source_basis/module_ao/parallel_orbitals.h"
 
 #include "source_estate/module_charge/charge.h" // mohan add 20251024
+#include "source_estate/module_dm/density_matrix.h" // mohan add 20251103
 
 namespace hsolver
 {
 
-template <typename T, typename Device = base_device::DEVICE_CPU>
+template <typename TK, typename Device = base_device::DEVICE_CPU>
 class HSolverLCAO
 {
   public:
     HSolverLCAO(const Parallel_Orbitals* ParaV_in, std::string method_in) : ParaV(ParaV_in), method(method_in) {};
 
-    void solve(hamilt::Hamilt<T>* pHamilt,
-               psi::Psi<T>& psi,
+    void solve(hamilt::Hamilt<TK>* pHamilt,
+               psi::Psi<TK>& psi,
                elecstate::ElecState* pes,
+			   elecstate::DensityMatrix<TK, double>& dm, // mohan add 2025-11-03
 			   Charge &chr, // charge density
 			   const int nspin,
 			   const bool skip_charge);
 
   private:
-    void hamiltSolvePsiK(hamilt::Hamilt<T>* hm, psi::Psi<T>& psi, double* eigenvalue); // for kpar_lcao == 1
+    void hamiltSolvePsiK(hamilt::Hamilt<TK>* hm, psi::Psi<TK>& psi, double* eigenvalue); // for kpar_lcao == 1
 
-    void parakSolve(hamilt::Hamilt<T>* pHamilt, psi::Psi<T>& psi, elecstate::ElecState* pes, int kpar); // for kpar_lcao > 1
+    void parakSolve(hamilt::Hamilt<TK>* pHamilt, psi::Psi<TK>& psi, elecstate::ElecState* pes, int kpar); // for kpar_lcao > 1
 
     // The solving algorithm using cusolver is different from others, so a separate function is needed
-    void parakSolve_cusolver(hamilt::Hamilt<T>* pHamilt,
-                             psi::Psi<T>& psi,
+    void parakSolve_cusolver(hamilt::Hamilt<TK>* pHamilt,
+                             psi::Psi<TK>& psi,
                              elecstate::ElecState* pes);
 
     const Parallel_Orbitals* ParaV;
diff --git a/source/source_io/cal_ldos.cpp b/source/source_io/cal_ldos.cpp
index f12fc9e391..009337c254 100644
--- a/source/source_io/cal_ldos.cpp
+++ b/source/source_io/cal_ldos.cpp
@@ -12,10 +12,16 @@ namespace ModuleIO
 
 #ifdef __LCAO
 template <typename T>
-void Cal_ldos<T>::cal_ldos_lcao(const elecstate::ElecStateLCAO<T>* pelec,
-                                const psi::Psi<T>& psi,
-                                const Parallel_Grid& pgrid,
-                                const UnitCell& ucell)
+void Cal_ldos<T>::cal_ldos_lcao(
+        const elecstate::Efermi &eferm, // mohan add 2025-11-02
+        const Charge &chr, // mohan add add 2025-11-02
+        const LCAO_domain::Setup_DM<T> &dmat, // mohan add 2025-11-02 
+		const K_Vectors &kv, // k points, mohan add 2025-11-02
+        const ModuleBase::matrix &ekb, // mohan add 2025-11-02
+        const ModuleBase::matrix &wg, // mohan add 2025-11-02
+		const psi::Psi<T>& psi,
+		const Parallel_Grid& pgrid,
+		const UnitCell& ucell)
 {
     for (int ie = 0; ie < PARAM.inp.stm_bias[2]; ie++)
     {
@@ -25,38 +31,38 @@ void Cal_ldos<T>::cal_ldos_lcao(const elecstate::ElecStateLCAO<T>* pelec,
         const double emax = en > 0 ? en : 0;
 
         // calculate weight (for bands not in the range, weight is zero)
-        ModuleBase::matrix weight(pelec->ekb.nr, pelec->ekb.nc);
-        for (int ik = 0; ik < pelec->ekb.nr; ++ik)
+        ModuleBase::matrix weight(ekb.nr, ekb.nc);
+        for (int ik = 0; ik < ekb.nr; ++ik)
         {
-            const double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]);
+            const double efermi = eferm.get_efval(kv.isk[ik]);
 
-            for (int ib = 0; ib < pelec->ekb.nc; ib++)
+            for (int ib = 0; ib < ekb.nc; ib++)
             {
-                const double eigenval = (pelec->ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV;
+                const double eigenval = (ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV;
                 if (eigenval >= emin && eigenval <= emax)
                 {
-                    weight(ik, ib) = en > 0 ? pelec->klist->wk[ik] - pelec->wg(ik, ib) : pelec->wg(ik, ib);
+                    weight(ik, ib) = en > 0 ? kv.wk[ik] - wg(ik, ib) : wg(ik, ib);
                 }
             }
         }
 
         // calculate dm-like for ldos
         const int nspin_dm = PARAM.inp.nspin == 2 ? 2 : 1;
-        elecstate::DensityMatrix<T, double> dm_ldos(pelec->DM->get_paraV_pointer(),
+        elecstate::DensityMatrix<T, double> dm_ldos(dmat.dm->get_paraV_pointer(),
                                                     nspin_dm,
-                                                    pelec->klist->kvec_d,
-                                                    pelec->klist->get_nks() / nspin_dm);
+                                                    kv.kvec_d,
+                                                    kv.get_nks() / nspin_dm);
 
-        elecstate::cal_dm_psi(pelec->DM->get_paraV_pointer(), weight, psi, dm_ldos);
-        dm_ldos.init_DMR(*(pelec->DM->get_DMR_pointer(1)));
+        elecstate::cal_dm_psi(dmat.dm->get_paraV_pointer(), weight, psi, dm_ldos);
+        dm_ldos.init_DMR(*(dmat.dm->get_DMR_pointer(1)));
         dm_ldos.cal_DMR();
 
         // allocate ldos space
-        std::vector<double> ldos_space(PARAM.inp.nspin * pelec->charge->nrxx);
+        std::vector<double> ldos_space(PARAM.inp.nspin * chr.nrxx);
         double** ldos = new double*[PARAM.inp.nspin];
         for (int is = 0; is < PARAM.inp.nspin; ++is)
         {
-            ldos[is] = &ldos_space[is * pelec->charge->nrxx];
+            ldos[is] = &ldos_space[is * chr.nrxx];
         }
 
     // calculate ldos
@@ -66,7 +72,7 @@ void Cal_ldos<T>::cal_ldos_lcao(const elecstate::ElecStateLCAO<T>* pelec,
         // ldos[0] += ldos[1] for nspin_dm == 2
         if (nspin_dm == 2)
         {
-            BlasConnector::axpy(pelec->charge->nrxx, 1.0, ldos[1], 1, ldos[0], 1);
+            BlasConnector::axpy(chr.nrxx, 1.0, ldos[1], 1, ldos[0], 1);
         }
 
         // write ldos to cube file
diff --git a/source/source_io/cal_ldos.h b/source/source_io/cal_ldos.h
index a18ac0f1d0..8a16046b76 100644
--- a/source/source_io/cal_ldos.h
+++ b/source/source_io/cal_ldos.h
@@ -4,6 +4,12 @@
 #include "source_estate/elecstate_lcao.h"
 #include "source_estate/elecstate_pw.h"
 
+#include "source_estate/fp_energy.h" // eferm
+#include "source_estate/module_charge/charge.h" // chr
+#include "source_lcao/setup_dm.h" // Setup_DM
+#include "source_cell/klist.h" // K_Vectors
+#include "source_base/matrix.h" // matrix
+
 namespace ModuleIO
 {
 template <typename T>
@@ -13,10 +19,17 @@ class Cal_ldos
     Cal_ldos(){};
     ~Cal_ldos(){};
 
-    static void cal_ldos_lcao(const elecstate::ElecStateLCAO<T>* pelec,
-                              const psi::Psi<T>& psi,
-                              const Parallel_Grid& pgrid,
-                              const UnitCell& ucell);
+    static void cal_ldos_lcao(
+        const elecstate::Efermi &eferm, // mohan add 2025-11-02
+        const Charge &chr, // mohan add add 2025-11-02
+        const LCAO_domain::Setup_DM<T> &dmat, // mohan add 2025-11-02 
+		const K_Vectors &kv, // k points, mohan add 2025-11-02
+        const ModuleBase::matrix &ekb, // mohan add 2025-11-02
+        const ModuleBase::matrix &wg, // mohan add 2025-11-02
+		const psi::Psi<T>& psi,
+		const Parallel_Grid& pgrid,
+		const UnitCell& ucell);
+
 }; // namespace Cal_ldos
 
 void cal_ldos_pw(const elecstate::ElecStatePW<std::complex<double>>* pelec,
@@ -72,4 +85,4 @@ void trilinear_interpolate(const std::vector<std::vector<int>>& points,
 
 } // namespace ModuleIO
 
-#endif // CAL_LDOS_H
\ No newline at end of file
+#endif // CAL_LDOS_H
diff --git a/source/source_io/ctrl_iter_lcao.cpp b/source/source_io/ctrl_iter_lcao.cpp
index d0b0c1bd0a..c772d162e6 100644
--- a/source/source_io/ctrl_iter_lcao.cpp
+++ b/source/source_io/ctrl_iter_lcao.cpp
@@ -13,7 +13,8 @@ template <typename TK, typename TR>
 void ctrl_iter_lcao(UnitCell& ucell, // unit cell *
         const Input_para& inp, // input parameters *
 		K_Vectors& kv, // k points *
-		elecstate::ElecStateLCAO<TK>* pelec, // electronic info * 
+		elecstate::ElecState* pelec, // electronic info * 
+        elecstate::DensityMatrix<TK, double>& dm, // density matrix, mohan add 2025-11-03
 		Parallel_Orbitals& pv, // parallel orbital info *
 		Grid_Driver& gd, // adjacent atom info *
 		psi::Psi<TK>* psi, // wave functions *
@@ -48,9 +49,9 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell *
         if (GlobalC::exx_info.info_global.cal_exx)
         {
             GlobalC::exx_info.info_ri.real_number ?
-              exx_nao.exd->exx_iter_finish(kv, ucell, *p_hamilt, *pelec,
+              exx_nao.exd->exx_iter_finish(kv, ucell, *p_hamilt, *pelec, &dm, 
                 *p_chgmix, scf_ene_thr, iter, istep, conv_esolver) :
-              exx_nao.exc->exx_iter_finish(kv, ucell, *p_hamilt, *pelec,
+              exx_nao.exc->exx_iter_finish(kv, ucell, *p_hamilt, *pelec, &dm,
                 *p_chgmix, scf_ene_thr, iter, istep, conv_esolver);
         }
     }
@@ -68,7 +69,7 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell *
 
             deepks_interface.out_deepks_labels(pelec->f_en.etot, kv.get_nks(),
               ucell.nat, PARAM.globalv.nlocal, pelec->ekb, kv.kvec_d,
-              ucell, orb, gd, &pv, *psi, pelec->get_DM(),
+              ucell, orb, gd, &pv, *psi, &dm,
               p_hamilt, iter, conv_esolver, GlobalV::MY_RANK, GlobalV::ofs_running);
         }
     }
@@ -81,7 +82,8 @@ void ctrl_iter_lcao(UnitCell& ucell, // unit cell *
 template void ctrl_iter_lcao<double, double>(UnitCell& ucell, // unit cell *
         const Input_para& inp, // input parameters *
 		K_Vectors& kv, // k points *
-		elecstate::ElecStateLCAO<double>* pelec, // electronic info * 
+		elecstate::ElecState* pelec, // electronic info * 
+        elecstate::DensityMatrix<double, double>& dm, // density matrix, mohan add 2025-11-03
 		Parallel_Orbitals& pv, // parallel orbital info *
 		Grid_Driver& gd, // adjacent atom info *
 		psi::Psi<double>* psi, // wave functions *
@@ -100,7 +102,8 @@ template void ctrl_iter_lcao<double, double>(UnitCell& ucell, // unit cell *
 template void ctrl_iter_lcao<std::complex<double>, double>(UnitCell& ucell, // unit cell *
         const Input_para& inp, // input parameters *
 		K_Vectors& kv, // k points *
-		elecstate::ElecStateLCAO<std::complex<double>>* pelec, // electronic info * 
+		elecstate::ElecState* pelec, // electronic info * 
+        elecstate::DensityMatrix<std::complex<double>, double>& dm, // density matrix, mohan add 2025-11-03
 		Parallel_Orbitals& pv, // parallel orbital info *
 		Grid_Driver& gd, // adjacent atom info *
 		psi::Psi<std::complex<double>>* psi, // wave functions *
@@ -119,7 +122,8 @@ template void ctrl_iter_lcao<std::complex<double>, double>(UnitCell& ucell, // u
 template void ctrl_iter_lcao<std::complex<double>, std::complex<double>>(UnitCell& ucell, // unit cell *
         const Input_para& inp, // input parameters *
 		K_Vectors& kv, // k points *
-		elecstate::ElecStateLCAO<std::complex<double>>* pelec, // electronic info * 
+		elecstate::ElecState* pelec, // electronic info * 
+        elecstate::DensityMatrix<std::complex<double>, double>& dm, // density matrix, mohan add 2025-11-03
 		Parallel_Orbitals& pv, // parallel orbital info *
 		Grid_Driver& gd, // adjacent atom info *
 		psi::Psi<std::complex<double>>* psi, // wave functions *
diff --git a/source/source_io/ctrl_iter_lcao.h b/source/source_io/ctrl_iter_lcao.h
index c94afac122..6297d4f1d0 100644
--- a/source/source_io/ctrl_iter_lcao.h
+++ b/source/source_io/ctrl_iter_lcao.h
@@ -18,7 +18,8 @@ template <typename TK, typename TR>
 void ctrl_iter_lcao(UnitCell& ucell, // unit cell *
         const Input_para& inp, // input parameters *
 		K_Vectors& kv, // k points *
-		elecstate::ElecStateLCAO<TK>* pelec, // electronic info * 
+		elecstate::ElecState* pelec, // electronic info * 
+        elecstate::DensityMatrix<TK, double>& dm, // density matrix, mohan add 2025-11-03
 		Parallel_Orbitals& pv, // parallel orbital info *
 		Grid_Driver& gd, // adjacent atom info *
 		psi::Psi<TK>* psi, // wave functions *
diff --git a/source/source_io/ctrl_runner_lcao.cpp b/source/source_io/ctrl_runner_lcao.cpp
index 9cf550830c..9a83e9f3d2 100644
--- a/source/source_io/ctrl_runner_lcao.cpp
+++ b/source/source_io/ctrl_runner_lcao.cpp
@@ -16,7 +16,8 @@ template <typename TK, typename TR>
 void ctrl_runner_lcao(UnitCell& ucell,      // unitcell
         const Input_para &inp,              // input
 		K_Vectors &kv,                      // k-point
-		elecstate::ElecStateLCAO<TK>* pelec,// electronic info
+		elecstate::ElecState* pelec,// electronic info
+        const LCAO_domain::Setup_DM<TK> &dmat, // mohan add 2025-11-02
 		Parallel_Orbitals &pv,              // orbital info
         Parallel_Grid &pgrid,               // grid info
 		Grid_Driver &gd,                    // search for adjacent atoms
@@ -44,7 +45,8 @@ void ctrl_runner_lcao(UnitCell& ucell,      // unitcell
 	// 2) out ldos
 	if (inp.out_ldos[0])
     {
-        ModuleIO::Cal_ldos<TK>::cal_ldos_lcao(pelec, psi[0], pgrid, ucell);
+        ModuleIO::Cal_ldos<TK>::cal_ldos_lcao(pelec->eferm, chr, dmat, kv, 
+          pelec->ekb, pelec->wg, psi[0], pgrid, ucell);
     }
 
     // 3) print out exchange-correlation potential
@@ -134,7 +136,8 @@ void ctrl_runner_lcao(UnitCell& ucell,      // unitcell
 template void ModuleIO::ctrl_runner_lcao<double, double>(UnitCell& ucell,      // unitcell
         const Input_para &inp,              // input
 		K_Vectors &kv,                      // k-point
-		elecstate::ElecStateLCAO<double>* pelec,// electronic info
+		elecstate::ElecState* pelec,// electronic info
+        const LCAO_domain::Setup_DM<double> &dmat, // mohan add 2025-11-02
 		Parallel_Orbitals &pv,              // orbital info
         Parallel_Grid &pgrid,               // grid info
 		Grid_Driver &gd,                    // search for adjacent atoms
@@ -154,7 +157,8 @@ template void ModuleIO::ctrl_runner_lcao<double, double>(UnitCell& ucell,      /
 template void ctrl_runner_lcao<std::complex<double>, double>(UnitCell& ucell,      // unitcell
         const Input_para &inp,              // input
 		K_Vectors &kv,                      // k-point
-		elecstate::ElecStateLCAO<std::complex<double>>* pelec,// electronic info
+		elecstate::ElecState* pelec,// electronic info
+        const LCAO_domain::Setup_DM<std::complex<double>> &dmat, // mohan add 2025-11-02
 		Parallel_Orbitals &pv,              // orbital info
         Parallel_Grid &pgrid,               // grid info
 		Grid_Driver &gd,                    // search for adjacent atoms
@@ -174,7 +178,8 @@ template void ctrl_runner_lcao<std::complex<double>, double>(UnitCell& ucell,
 template void ctrl_runner_lcao<std::complex<double>, std::complex<double>>(UnitCell& ucell,      // unitcell
         const Input_para &inp,              // input
 		K_Vectors &kv,                      // k-point
-		elecstate::ElecStateLCAO<std::complex<double>>* pelec,// electronic info
+		elecstate::ElecState* pelec,// electronic info
+        const LCAO_domain::Setup_DM<std::complex<double>> &dmat, // mohan add 2025-11-02
 		Parallel_Orbitals &pv,              // orbital info
         Parallel_Grid &pgrid,               // grid info
 		Grid_Driver &gd,                    // search for adjacent atoms
diff --git a/source/source_io/ctrl_runner_lcao.h b/source/source_io/ctrl_runner_lcao.h
index 2b57c1800f..95c4fe4658 100644
--- a/source/source_io/ctrl_runner_lcao.h
+++ b/source/source_io/ctrl_runner_lcao.h
@@ -3,11 +3,12 @@
 
 #include "source_cell/unitcell.h" // use UnitCell
 #include "source_cell/klist.h" // use K_Vectors
-#include "source_estate/elecstate_lcao.h" // use elecstate::ElecStateLCAO<TK> 
+#include "source_estate/elecstate.h" // use elecstate::ElecStateLCAO<TK> 
 #include "source_psi/psi.h" // use Psi<TK>
 #include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO<TK, TR>
 #include "source_basis/module_nao/two_center_bundle.h" // use TwoCenterBundle
 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251018
+#include "source_lcao/setup_dm.h" // for density matrix, mohan add 20251103
 
 namespace ModuleIO
 {
@@ -16,7 +17,8 @@ template <typename TK, typename TR>
 void ctrl_runner_lcao(UnitCell& ucell,      // unitcell
         const Input_para &inp,              // input
 		K_Vectors &kv,                      // k-point
-		elecstate::ElecStateLCAO<TK>* pelec,// electronic info
+		elecstate::ElecState* pelec,// electronic info
+        const LCAO_domain::Setup_DM<TK> &dmat, // mohan add 2025-11-02
 		Parallel_Orbitals &pv,              // orbital info
         Parallel_Grid &pgrid,               // grid info
 		Grid_Driver &gd,                    // search for adjacent atoms
diff --git a/source/source_io/ctrl_scf_lcao.cpp b/source/source_io/ctrl_scf_lcao.cpp
index c859cf3f74..f66fd2720b 100644
--- a/source/source_io/ctrl_scf_lcao.cpp
+++ b/source/source_io/ctrl_scf_lcao.cpp
@@ -36,7 +36,8 @@ template <typename TK, typename TR>
 void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
         const Input_para& inp,
 		K_Vectors& kv,
-		elecstate::ElecStateLCAO<TK>* pelec, 
+		elecstate::ElecState* pelec, 
+		elecstate::DensityMatrix<TK,double>* dm, // mohan add 2025-11-04
 		Parallel_Orbitals& pv,
 		Grid_Driver& gd,
 		psi::Psi<TK>* psi,
@@ -115,11 +116,9 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
 	//------------------------------------------------------------------
     if(inp.out_dmr[0])
 	{
-		const auto& dmr_vector = pelec->get_DM()->get_DMR_vector();
-
         const int precision = inp.out_dmr[1];
 
-		ModuleIO::write_dmr(dmr_vector, precision, pv, out_app_flag,
+		ModuleIO::write_dmr(dm->get_DMR_vector(), precision, pv, out_app_flag,
 				ucell.get_iat2iwt(), ucell.nat, istep);
 	}
 
@@ -135,7 +134,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
 		}
 		const int precision = inp.out_dmk[1];
 
-		ModuleIO::write_dmk(pelec->get_DM()->get_DMK_vector(),
+		ModuleIO::write_dmk(dm->get_DMK_vector(),
 				precision, efermis, &(ucell), pv, istep);
 	}
 
@@ -196,7 +195,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
 			gd,
 			&pv,
 			*psi,
-			pelec->get_DM(),
+			dm,
 			p_ham_deepks,
             -1, // -1 when called in after scf
             true, // no used when after scf
@@ -300,7 +299,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
         ModuleIO::cal_mag(&pv,
                 p_hamilt,
                 kv,
-                pelec,
+                dm, // mohan add 2025-11-04
                 two_center_bundle,
                 orb,
                 ucell,
@@ -404,7 +403,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
     if (inp.rpa)
     {
         RPA_LRI<TK, double> rpa_lri_double(GlobalC::exx_info.info_ri);
-        rpa_lri_double.cal_postSCF_exx(*dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(pelec)->get_DM(),
+        rpa_lri_double.cal_postSCF_exx(*dm,
                                        MPI_COMM_WORLD,
                                        ucell,
                                        kv,
@@ -469,7 +468,8 @@ template void ModuleIO::ctrl_scf_lcao<double, double>(
         UnitCell& ucell, 
         const Input_para& inp,
 		K_Vectors& kv,
-		elecstate::ElecStateLCAO<double>* pelec, 
+		elecstate::ElecState* pelec, 
+        elecstate::DensityMatrix<double,double>* dm, // mohan add 2025-11-04
 		Parallel_Orbitals& pv,
 		Grid_Driver& gd,
 		psi::Psi<double>* psi,
@@ -492,7 +492,8 @@ template void ModuleIO::ctrl_scf_lcao<std::complex<double>, double>(
         UnitCell& ucell, 
         const Input_para& inp,
 		K_Vectors& kv,
-		elecstate::ElecStateLCAO<std::complex<double>>* pelec, 
+		elecstate::ElecState* pelec, 
+        elecstate::DensityMatrix<std::complex<double>,double>* dm, // mohan add 2025-11-04
 		Parallel_Orbitals& pv,
 		Grid_Driver& gd,
 		psi::Psi<std::complex<double>>* psi,
@@ -514,7 +515,8 @@ template void ModuleIO::ctrl_scf_lcao<std::complex<double>, std::complex<double>
         UnitCell& ucell, 
         const Input_para& inp,
 		K_Vectors& kv,
-		elecstate::ElecStateLCAO<std::complex<double>>* pelec, 
+		elecstate::ElecState* pelec, 
+        elecstate::DensityMatrix<std::complex<double>,double>* dm, // mohan add 2025-11-04
 		Parallel_Orbitals& pv,
 		Grid_Driver& gd,
 		psi::Psi<std::complex<double>>* psi,
diff --git a/source/source_io/ctrl_scf_lcao.h b/source/source_io/ctrl_scf_lcao.h
index 98ee5e18a5..9dd381cd2a 100644
--- a/source/source_io/ctrl_scf_lcao.h
+++ b/source/source_io/ctrl_scf_lcao.h
@@ -5,7 +5,7 @@
 
 #include "source_cell/unitcell.h" // use UnitCell
 #include "source_cell/klist.h" // use K_Vectors
-#include "source_estate/elecstate_lcao.h" // use elecstate::ElecStateLCAO<TK> 
+#include "source_estate/elecstate.h" // use elecstate::ElecStateLCAO<TK> 
 #include "source_psi/psi.h" // use Psi<TK>
 #include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO<TK, TR>
 #include "source_basis/module_nao/two_center_bundle.h" // use TwoCenterBundle
@@ -15,6 +15,7 @@
 
 #include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251008
 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008
+#include "source_estate/module_dm/density_matrix.h" // mohan add 2025-11-04
 
 namespace ModuleIO
 {
@@ -23,7 +24,8 @@ namespace ModuleIO
 		void ctrl_scf_lcao(UnitCell& ucell, 
 				const Input_para& inp,
 				K_Vectors& kv,
-				elecstate::ElecStateLCAO<TK>* pelec, 
+				elecstate::ElecState* pelec, 
+                elecstate::DensityMatrix<TK,double> *dm, // mohan add 2025-11-04
 				Parallel_Orbitals& pv,
 				Grid_Driver& gd,
 				psi::Psi<TK>* psi,
diff --git a/source/source_io/get_wf_lcao.cpp b/source/source_io/get_wf_lcao.cpp
index 7d6cd0d15c..6cafa72d60 100644
--- a/source/source_io/get_wf_lcao.cpp
+++ b/source/source_io/get_wf_lcao.cpp
@@ -480,6 +480,7 @@ template int Get_wf_lcao::set_wfc_grid(const int naroc[2],
                                        const double* in,
                                        double** out,
                                        const std::vector<int>& trace_lo);
+
 template int Get_wf_lcao::set_wfc_grid(const int naroc[2],
                                        const int nb,
                                        const int dim0,
diff --git a/source/source_io/output_mulliken.h b/source/source_io/output_mulliken.h
index e20cea9e04..42c3cdc550 100644
--- a/source/source_io/output_mulliken.h
+++ b/source/source_io/output_mulliken.h
@@ -9,6 +9,7 @@
 #include "source_io/output_sk.h"
 #include "source_base/formatter.h"
 #include "source_lcao/module_operator_lcao/dspin_lcao.h"
+#include "source_estate/module_dm/density_matrix.h" // mohan add 2025-11-04
 
 #include <map>
 #include <vector>
@@ -89,7 +90,7 @@ template <typename TK>
 void cal_mag(Parallel_Orbitals* pv,
              hamilt::Hamilt<TK>* p_ham,
              K_Vectors& kv,
-             elecstate::ElecState* pelec,
+             elecstate::DensityMatrix<TK,double>* dm, // mohan add 2025-11-04
              const TwoCenterBundle& two_center_bundle,
              const LCAO_Orbitals& orb,
              UnitCell& ucell,
@@ -103,10 +104,7 @@ void cal_mag(Parallel_Orbitals* pv,
         auto cell_index
             = CellIndex(ucell.get_atomLabels(), ucell.get_atomCounts(), ucell.get_lnchiCounts(), PARAM.inp.nspin);
         auto out_s_k = ModuleIO::Output_Sk<TK>(p_ham, pv, PARAM.inp.nspin, kv.get_nks());
-        auto out_dm_k = ModuleIO::Output_DMK<TK>(dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(pelec)->get_DM(),
-                                                pv,
-                                                PARAM.inp.nspin,
-                                                kv.get_nks());
+        auto out_dm_k = ModuleIO::Output_DMK<TK>(dm, pv, PARAM.inp.nspin, kv.get_nks());
 
         auto mulp = ModuleIO::Output_Mulliken<TK>(&(out_s_k), &(out_dm_k), pv, &cell_index, kv.isk, PARAM.inp.nspin);
         auto atom_chg = mulp.get_atom_chg();
@@ -127,8 +125,7 @@ void cal_mag(Parallel_Orbitals* pv,
     {
         std::vector<std::vector<double>> atom_mag(ucell.nat, std::vector<double>(PARAM.inp.nspin, 0.0));
         std::vector<ModuleBase::Vector3<int>> constrain(ucell.nat, ModuleBase::Vector3<int>(1, 1, 1));
-        const hamilt::HContainer<double>* dmr
-            = dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(pelec)->get_DM()->get_DMR_pointer(1);
+        const hamilt::HContainer<double>* dmr = dm->get_DMR_pointer(1);
         std::vector<double> moments;
         std::vector<double> mag_x(ucell.nat, 0.0);
         std::vector<double> mag_y(ucell.nat, 0.0);
@@ -144,9 +141,9 @@ void cal_mag(Parallel_Orbitals* pv,
 						&gd,
 						two_center_bundle.overlap_orb_onsite.get(),
 						orb.cutoffs());
-			dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(pelec)->get_DM()->switch_dmr(2);
+			dm->switch_dmr(2);
 			moments = sc_lambda->cal_moment(dmr, constrain);
-			dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(pelec)->get_DM()->switch_dmr(0);
+			dm->switch_dmr(0);
 			delete sc_lambda;
 			//const std::vector<std::string> title = {"Total Magnetism (uB)", ""};
             //const std::vector<std::string> fmts = {"%-26s", "%20.10f"};
diff --git a/source/source_lcao/CMakeLists.txt b/source/source_lcao/CMakeLists.txt
index dbd2113bfb..5199625c41 100644
--- a/source/source_lcao/CMakeLists.txt
+++ b/source/source_lcao/CMakeLists.txt
@@ -42,6 +42,7 @@ if(ENABLE_LCAO)
         LCAO_init_basis.cpp
         setup_exx.cpp
         setup_deepks.cpp
+        setup_dm.cpp
         rho_tau_lcao.cpp
 		record_adj.cpp
         center2_orb.cpp
diff --git a/source/source_lcao/FORCE.h b/source/source_lcao/FORCE.h
index 5eba250181..a80dd02f25 100644
--- a/source/source_lcao/FORCE.h
+++ b/source/source_lcao/FORCE.h
@@ -37,8 +37,9 @@ class Force_LCAO
                 const UnitCell& ucell,
                 const Grid_Driver& gd,
                 const psi::Psi<T>* psi,
-                const elecstate::ElecState* pelec,
-                ModuleBase::matrix& foverlap,
+				const elecstate::ElecState* pelec,
+				const elecstate::DensityMatrix<T, double>* dm, // mohan add 2025-11-04
+				ModuleBase::matrix& foverlap,
                 ModuleBase::matrix& ftvnl_dphi,
                 ModuleBase::matrix& fvnl_dbeta,
                 ModuleBase::matrix& fvl_dphi,
diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp
index af51e04ab2..54fe92ed9a 100644
--- a/source/source_lcao/FORCE_STRESS.cpp
+++ b/source/source_lcao/FORCE_STRESS.cpp
@@ -22,6 +22,34 @@
 #include "source_lcao/module_operator_lcao/dspin_lcao.h"
 #include "source_lcao/module_operator_lcao/nonlocal_new.h"
 
+
+// mohan add 2025-11-04
+template <>
+void assign_dmk_ptr<double>(
+    elecstate::DensityMatrix<double,double>* dm,
+    std::vector<std::vector<double>>*& dmk_d,
+    std::vector<std::vector<std::complex<double>>>*& dmk_c,
+    bool gamma_only_local
+) {
+    auto& dmk_tmp = dm->get_DMK_vector();
+    dmk_d = &dmk_tmp;
+    dmk_c = nullptr;
+}
+
+template <>
+void assign_dmk_ptr<std::complex<double>>(
+    elecstate::DensityMatrix<std::complex<double>,double>* dm,
+    std::vector<std::vector<double>>*& dmk_d,
+    std::vector<std::vector<std::complex<double>>>*& dmk_c,
+    bool gamma_only_local
+) {
+    auto& dmk_tmp = dm->get_DMK_vector();
+    dmk_c = &dmk_tmp;
+    dmk_d = nullptr;
+}
+
+
+
 template <typename T>
 Force_Stress_LCAO<T>::Force_Stress_LCAO(Record_adj& ra, const int nat_in) : RA(&ra), nat(nat_in)
 {
@@ -39,6 +67,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
                                           const Grid_Driver& gd,
                                           Parallel_Orbitals& pv,
                                           const elecstate::ElecState* pelec,
+                                          LCAO_domain::Setup_DM<T> &dmat, // mohan add 2025-11-03
                                           const psi::Psi<T>* psi,
                                           const TwoCenterBundle& two_center_bundle,
                                           const LCAO_Orbitals& orb,
@@ -133,7 +162,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
 
     //! atomic forces from integration (4 terms)
     this->integral_part(PARAM.globalv.gamma_only_local, isforce, isstress,
-                        ucell, gd, fsr, pelec, psi, foverlap, ftvnl_dphi,
+                        ucell, gd, fsr, pelec, dmat.dm, psi, foverlap, ftvnl_dphi, // add dmat.dm, mohan 20251104
                         fvnl_dbeta, fvl_dphi, soverlap, stvnl_dphi, svnl_dbeta,
                         svl_dphi, fvnl_dalpha, svnl_dalpha, deepks, 
                         two_center_bundle, orb, pv, kv);
@@ -144,18 +173,16 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
         hamilt::NonlocalNew<hamilt::OperatorLCAO<T, double>> tmp_nonlocal(nullptr,
           kv.kvec_d, nullptr, &ucell, orb.cutoffs(), &gd, two_center_bundle.overlap_orb_beta.get());
 
-        const auto* dm_p = dynamic_cast<const elecstate::ElecStateLCAO<T>*>(pelec)->get_DM();
-
         if (PARAM.inp.nspin == 2)
         {
-            const_cast<elecstate::DensityMatrix<T, double>*>(dm_p)->switch_dmr(1);
+            dmat.dm->switch_dmr(1);
         }
 
-        const hamilt::HContainer<double>* dmr = dm_p->get_DMR_pointer(1);
+        const hamilt::HContainer<double>* dmr = dmat.dm->get_DMR_pointer(1);
         tmp_nonlocal.cal_force_stress(isforce, isstress, dmr, fvnl_dbeta, svnl_dbeta);
         if (PARAM.inp.nspin == 2)
         {
-            const_cast<elecstate::DensityMatrix<T, double>*>(dm_p)->switch_dmr(0);
+            dmat.dm->switch_dmr(0);
         }
     }
     else if (PARAM.inp.nspin == 4)
@@ -166,12 +193,11 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
 
         // calculate temporary complex DMR for nonlocal force&stress
         // In fact, only SOC part need the imaginary part of DMR for correct force&stress
-        const auto* dm_p = dynamic_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(pelec)->get_DM();
-        hamilt::HContainer<std::complex<double>> tmp_dmr(dm_p->get_DMR_pointer(1)->get_paraV());
-        std::vector<int> ijrs = dm_p->get_DMR_pointer(1)->get_ijr_info();
+        hamilt::HContainer<std::complex<double>> tmp_dmr(dmat.dm->get_DMR_pointer(1)->get_paraV());
+        std::vector<int> ijrs = dmat.dm->get_DMR_pointer(1)->get_ijr_info();
         tmp_dmr.insert_ijrs(&ijrs);
         tmp_dmr.allocate();
-        dm_p->cal_DMR_full(&tmp_dmr);
+        dmat.dm->cal_DMR_full(&tmp_dmr);
         tmp_nonlocal.cal_force_stress(isforce, isstress, &tmp_dmr, fvnl_dbeta, svnl_dbeta);
     }
 
@@ -248,7 +274,13 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
         }
         if (PARAM.inp.dft_plus_u == 2)
         {
-            GlobalC::dftu.force_stress(ucell, gd, pelec, pv, fsr, force_dftu, stress_dftu, kv);
+			// GlobalC::dftu.force_stress(ucell, gd, pelec, pv, fsr, force_dftu, stress_dftu, kv);
+			// mohan modify 2025-11-03
+            std::vector<std::vector<double>>* dmk_d = nullptr;
+            std::vector<std::vector<std::complex<double>>>* dmk_c = nullptr;
+            // add a new template function
+            assign_dmk_ptr<T>(dmat.dm, dmk_d, dmk_c, PARAM.globalv.gamma_only_local);
+		    GlobalC::dftu.force_stress(ucell, gd, dmk_d, dmk_c, pv, fsr, force_dftu, stress_dftu, kv);
         }
         else
         {
@@ -287,16 +319,15 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
                                                                      two_center_bundle.overlap_orb_onsite.get(),
                                                                      orb.cutoffs());
 
-        const auto* dm_p = dynamic_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(pelec)->get_DM();
         if (PARAM.inp.nspin == 2)
         {
-            const_cast<elecstate::DensityMatrix<std::complex<double>, double>*>(dm_p)->switch_dmr(2);
+            dmat.dm->switch_dmr(2);
         }
-        const hamilt::HContainer<double>* dmr = dm_p->get_DMR_pointer(1);
+        const hamilt::HContainer<double>* dmr = dmat.dm->get_DMR_pointer(1);
         tmp_dspin.cal_force_stress(isforce, isstress, dmr, force_dspin, stress_dspin);
         if (PARAM.inp.nspin == 2)
         {
-            const_cast<elecstate::DensityMatrix<std::complex<double>, double>*>(dm_p)->switch_dmr(0);
+            dmat.dm->switch_dmr(0);
         }
     }
 
@@ -712,31 +743,32 @@ void Force_Stress_LCAO<T>::calForcePwPart(UnitCell& ucell,
 // overlap, kinetic, nonlocal pseudopotential, Local potential terms in force and stress
 template <>
 void Force_Stress_LCAO<double>::integral_part(const bool isGammaOnly,
-                                              const bool isforce,
-                                              const bool isstress,
-                                              const UnitCell& ucell,
-                                              const Grid_Driver& gd,
-                                              ForceStressArrays& fsr, // mohan add 2024-06-15
-                                              const elecstate::ElecState* pelec,
-                                              const psi::Psi<double>* psi,
-                                              ModuleBase::matrix& foverlap,
-                                              ModuleBase::matrix& ftvnl_dphi,
-                                              ModuleBase::matrix& fvnl_dbeta,
-                                              ModuleBase::matrix& fvl_dphi,
-                                              ModuleBase::matrix& soverlap,
-                                              ModuleBase::matrix& stvnl_dphi,
-                                              ModuleBase::matrix& svnl_dbeta,
-                                              ModuleBase::matrix& svl_dphi,
-                                              ModuleBase::matrix& fvnl_dalpha,
-                                              ModuleBase::matrix& svnl_dalpha,
-                                              Setup_DeePKS<double>& deepks,
-                                              const TwoCenterBundle& two_center_bundle,
-                                              const LCAO_Orbitals& orb,
-                                              const Parallel_Orbitals& pv,
-                                              const K_Vectors& kv)
+		const bool isforce,
+		const bool isstress,
+		const UnitCell& ucell,
+		const Grid_Driver& gd,
+		ForceStressArrays& fsr, // mohan add 2024-06-15
+		const elecstate::ElecState* pelec,
+		const elecstate::DensityMatrix<double, double>* dm, // mohan add 2025-11-04
+		const psi::Psi<double>* psi,
+		ModuleBase::matrix& foverlap,
+		ModuleBase::matrix& ftvnl_dphi,
+		ModuleBase::matrix& fvnl_dbeta,
+		ModuleBase::matrix& fvl_dphi,
+		ModuleBase::matrix& soverlap,
+		ModuleBase::matrix& stvnl_dphi,
+		ModuleBase::matrix& svnl_dbeta,
+		ModuleBase::matrix& svl_dphi,
+		ModuleBase::matrix& fvnl_dalpha,
+		ModuleBase::matrix& svnl_dalpha,
+		Setup_DeePKS<double>& deepks,
+		const TwoCenterBundle& two_center_bundle,
+		const LCAO_Orbitals& orb,
+		const Parallel_Orbitals& pv,
+		const K_Vectors& kv)
 {
 
-    flk.ftable(isforce, isstress, fsr, ucell, gd, psi, pelec,
+    flk.ftable(isforce, isstress, fsr, ucell, gd, psi, pelec, dm,
                foverlap, ftvnl_dphi, fvnl_dbeta, fvl_dphi,
                soverlap, stvnl_dphi, svnl_dbeta, svl_dphi,
                fvnl_dalpha, svnl_dalpha, deepks, two_center_bundle, orb, pv);
@@ -745,30 +777,31 @@ void Force_Stress_LCAO<double>::integral_part(const bool isGammaOnly,
 
 template <>
 void Force_Stress_LCAO<std::complex<double>>::integral_part(const bool isGammaOnly,
-                                                            const bool isforce,
-                                                            const bool isstress,
-                                                            const UnitCell& ucell,
-                                                            const Grid_Driver& gd,
-                                                            ForceStressArrays& fsr, // mohan add 2024-06-15
-                                                            const elecstate::ElecState* pelec,
-                                                            const psi::Psi<std::complex<double>>* psi,
-                                                            ModuleBase::matrix& foverlap,
-                                                            ModuleBase::matrix& ftvnl_dphi,
-                                                            ModuleBase::matrix& fvnl_dbeta,
-                                                            ModuleBase::matrix& fvl_dphi,
-                                                            ModuleBase::matrix& soverlap,
-                                                            ModuleBase::matrix& stvnl_dphi,
-                                                            ModuleBase::matrix& svnl_dbeta,
-                                                            ModuleBase::matrix& svl_dphi,
-                                                            ModuleBase::matrix& fvnl_dalpha,
-                                                            ModuleBase::matrix& svnl_dalpha,
-                                                            Setup_DeePKS<std::complex<double>>& deepks,
-                                                            const TwoCenterBundle& two_center_bundle,
-                                                            const LCAO_Orbitals& orb,
-                                                            const Parallel_Orbitals& pv,
-                                                            const K_Vectors& kv)
+		const bool isforce,
+		const bool isstress,
+		const UnitCell& ucell,
+		const Grid_Driver& gd,
+		ForceStressArrays& fsr, // mohan add 2024-06-15
+		const elecstate::ElecState* pelec,
+		const elecstate::DensityMatrix<std::complex<double>, double>* dm, // mohan add 2025-11-04
+		const psi::Psi<std::complex<double>>* psi,
+		ModuleBase::matrix& foverlap,
+		ModuleBase::matrix& ftvnl_dphi,
+		ModuleBase::matrix& fvnl_dbeta,
+		ModuleBase::matrix& fvl_dphi,
+		ModuleBase::matrix& soverlap,
+		ModuleBase::matrix& stvnl_dphi,
+		ModuleBase::matrix& svnl_dbeta,
+		ModuleBase::matrix& svl_dphi,
+		ModuleBase::matrix& fvnl_dalpha,
+		ModuleBase::matrix& svnl_dalpha,
+		Setup_DeePKS<std::complex<double>>& deepks,
+		const TwoCenterBundle& two_center_bundle,
+		const LCAO_Orbitals& orb,
+		const Parallel_Orbitals& pv,
+		const K_Vectors& kv)
 {
-    flk.ftable(isforce, isstress, fsr, ucell, gd, psi, pelec,
+    flk.ftable(isforce, isstress, fsr, ucell, gd, psi, pelec, dm,
                foverlap, ftvnl_dphi, fvnl_dbeta, fvl_dphi,
                soverlap, stvnl_dphi, svnl_dbeta, svl_dphi,
                fvnl_dalpha, svnl_dalpha, deepks,
diff --git a/source/source_lcao/FORCE_STRESS.h b/source/source_lcao/FORCE_STRESS.h
index dd7d464b94..b62dd1460b 100644
--- a/source/source_lcao/FORCE_STRESS.h
+++ b/source/source_lcao/FORCE_STRESS.h
@@ -16,6 +16,7 @@
 #include "force_stress_arrays.h"
 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008
 #include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251010
+#include "source_lcao/setup_dm.h" // mohan add 2025-11-03
 
 
 template <typename T>
@@ -38,6 +39,7 @@ class Force_Stress_LCAO
                         const Grid_Driver& gd,
                         Parallel_Orbitals& pv,
                         const elecstate::ElecState* pelec,
+                        LCAO_domain::Setup_DM<T> &dmat, // mohan add 2025-11-03
                         const psi::Psi<T>* psi,
                         const TwoCenterBundle& two_center_bundle,
                         const LCAO_Orbitals& orb,
@@ -79,8 +81,9 @@ class Force_Stress_LCAO
                        const UnitCell& ucell,
                        const Grid_Driver& gd,
                        ForceStressArrays& fsr, // mohan add 2024-06-15
-                       const elecstate::ElecState* pelec,
-                       const psi::Psi<T>* psi,
+					   const elecstate::ElecState* pelec,
+					   const elecstate::DensityMatrix<T, double>* dm, // mohan add 2025-11-04
+					   const psi::Psi<T>* psi,
                        ModuleBase::matrix& foverlap,
                        ModuleBase::matrix& ftvnl_dphi,
                        ModuleBase::matrix& fvnl_dbeta,
@@ -115,4 +118,13 @@ class Force_Stress_LCAO
 template <typename T>
 double Force_Stress_LCAO<T>::force_invalid_threshold_ev = 0.00;
 
+// only for DFT+U, mohan add 2025-11-04
+template <typename T>
+void assign_dmk_ptr(
+    elecstate::DensityMatrix<T,double>* dm,
+    std::vector<std::vector<double>>*& dmk_d,
+    std::vector<std::vector<std::complex<double>>>*& dmk_c,
+    bool gamma_only_local
+);
+
 #endif
diff --git a/source/source_lcao/FORCE_gamma.cpp b/source/source_lcao/FORCE_gamma.cpp
index 763a3e3321..fb9249c4b0 100644
--- a/source/source_lcao/FORCE_gamma.cpp
+++ b/source/source_lcao/FORCE_gamma.cpp
@@ -158,9 +158,10 @@ void Force_LCAO<double>::ftable(const bool isforce,
                                 const UnitCell& ucell,
                                 const Grid_Driver& gd,
                                 const psi::Psi<double>* psi,
-                                const elecstate::ElecState* pelec,
-                                ModuleBase::matrix& foverlap,
-                                ModuleBase::matrix& ftvnl_dphi,
+								const elecstate::ElecState* pelec,
+								const elecstate::DensityMatrix<double, double>* dm, // mohan add 2025-11-04
+								ModuleBase::matrix& foverlap,
+								ModuleBase::matrix& ftvnl_dphi,
                                 ModuleBase::matrix& fvnl_dbeta,
                                 ModuleBase::matrix& fvl_dphi,
                                 ModuleBase::matrix& soverlap,
@@ -179,10 +180,6 @@ void Force_LCAO<double>::ftable(const bool isforce,
     ModuleBase::TITLE("Forces", "ftable");
     ModuleBase::timer::tick("Forces", "ftable");
 
-    // get DM
-    const elecstate::DensityMatrix<double, double>* dm
-        = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec)->get_DM();
-
     this->ParaV = dm->get_paraV_pointer();
 
     // allocate DSloc_x, DSloc_y, DSloc_z
diff --git a/source/source_lcao/FORCE_k.cpp b/source/source_lcao/FORCE_k.cpp
index 1e19a53a20..7a774481b7 100644
--- a/source/source_lcao/FORCE_k.cpp
+++ b/source/source_lcao/FORCE_k.cpp
@@ -190,35 +190,33 @@ void Force_LCAO<std::complex<double>>::finish_ftable(ForceStressArrays& fsr)
 // be called in Force_LCAO::start_force_calculation
 template <>
 void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
-                                              const bool isstress,
-                                              ForceStressArrays& fsr, // mohan add 2024-06-15
-                                              const UnitCell& ucell,
-                                              const Grid_Driver& gd,
-                                              const psi::Psi<std::complex<double>>* psi,
-                                              const elecstate::ElecState* pelec,
-                                              ModuleBase::matrix& foverlap,
-                                              ModuleBase::matrix& ftvnl_dphi,
-                                              ModuleBase::matrix& fvnl_dbeta,
-                                              ModuleBase::matrix& fvl_dphi,
-                                              ModuleBase::matrix& soverlap,
-                                              ModuleBase::matrix& stvnl_dphi,
-                                              ModuleBase::matrix& svnl_dbeta,
-                                              ModuleBase::matrix& svl_dphi,
-                                              ModuleBase::matrix& fvnl_dalpha,
-                                              ModuleBase::matrix& svnl_dalpha,
-                                              Setup_DeePKS<std::complex<double>>& deepks,
-                                              const TwoCenterBundle& two_center_bundle,
-                                              const LCAO_Orbitals& orb,
-                                              const Parallel_Orbitals& pv,
-                                              const K_Vectors* kv,
-                                              Record_adj* ra)
+		const bool isstress,
+		ForceStressArrays& fsr, // mohan add 2024-06-15
+		const UnitCell& ucell,
+		const Grid_Driver& gd,
+		const psi::Psi<std::complex<double>>* psi,
+		const elecstate::ElecState* pelec,
+		const elecstate::DensityMatrix<std::complex<double>, double>* dm, // mohan add 2025-11-04
+		ModuleBase::matrix& foverlap,
+		ModuleBase::matrix& ftvnl_dphi,
+		ModuleBase::matrix& fvnl_dbeta,
+		ModuleBase::matrix& fvl_dphi,
+		ModuleBase::matrix& soverlap,
+		ModuleBase::matrix& stvnl_dphi,
+		ModuleBase::matrix& svnl_dbeta,
+		ModuleBase::matrix& svl_dphi,
+		ModuleBase::matrix& fvnl_dalpha,
+		ModuleBase::matrix& svnl_dalpha,
+		Setup_DeePKS<std::complex<double>>& deepks,
+		const TwoCenterBundle& two_center_bundle,
+		const LCAO_Orbitals& orb,
+		const Parallel_Orbitals& pv,
+		const K_Vectors* kv,
+		Record_adj* ra)
 {
     ModuleBase::TITLE("Forces", "ftable");
     ModuleBase::timer::tick("Forces", "ftable");
 
-    elecstate::DensityMatrix <std::complex<double>, double>* dm
-        = dynamic_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(pelec)->get_DM();
-
     this->allocate(ucell,
                    gd,
                    pv,
diff --git a/source/source_lcao/module_deltaspin/cal_mw.cpp b/source/source_lcao/module_deltaspin/cal_mw.cpp
index 7f6a49e207..becdd82e7c 100644
--- a/source/source_lcao/module_deltaspin/cal_mw.cpp
+++ b/source/source_lcao/module_deltaspin/cal_mw.cpp
@@ -19,14 +19,16 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mi_lcao(const int&
     ModuleBase::timer::tick("spinconstrain::SpinConstrain", "cal_mi_lcao");
     // calculate MW from lambda in real space projection method
     this->zero_Mi();
-    const hamilt::HContainer<double>* dmr
-        = static_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->get_DMR_pointer(1);
+    const hamilt::HContainer<double>* dmr = this->dm_->get_DMR_pointer(1);
     std::vector<double> moments;
     if(PARAM.inp.nspin==2)
     {
-        static_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->switch_dmr(2);
+        this->dm_->switch_dmr(2);
+
         moments = static_cast<hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, double>>*>(this->p_operator)->cal_moment(dmr, this->get_constrain());
-        static_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->switch_dmr(0);
+
+        this->dm_->switch_dmr(0);
+
         for(int iat=0;iat<this->Mi_.size();iat++)
         {
             this->Mi_[iat].x = 0.0;
@@ -167,4 +169,4 @@ void spinconstrain::SpinConstrain<double>::set_operator(
     hamilt::Operator<double>* op_in)
 {
     this->p_operator = op_in;
-}
\ No newline at end of file
+}
diff --git a/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp b/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp
index 515a89ab9c..e8bb81db23 100644
--- a/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp
+++ b/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp
@@ -128,7 +128,9 @@ void spinconstrain::SpinConstrain<std::complex<double>>::calculate_delta_hcc(std
 }
 
 template <>
-void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(int i_step, const ModuleBase::Vector3<double>* delta_lambda)
+void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
+		int i_step, 
+		const ModuleBase::Vector3<double>* delta_lambda)
 {
     ModuleBase::TITLE("spinconstrain::SpinConstrain", "cal_mw_from_lambda");
     ModuleBase::timer::tick("spinconstrain::SpinConstrain", "cal_mw_from_lambda");
@@ -152,7 +154,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(int
         }
         // diagonalization without update charge
         // mohan add two parameters charge and nspin, 2025-10-24
-        hsolver_t.solve(hamilt_t, psi_t[0], this->pelec, *this->pelec->charge, PARAM.inp.nspin, true);
+        hsolver_t.solve(hamilt_t, psi_t[0], this->pelec, *this->dm_, *this->pelec->charge, PARAM.inp.nspin, true);
         elecstate::calculate_weights(this->pelec->ekb,
                                      this->pelec->wg,
                                      this->pelec->klist,
@@ -161,10 +163,11 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(int
                                      this->pelec->nelec_spin,
                                      this->pelec->skip_weights);
         elecstate::calEBand(this->pelec->ekb,this->pelec->wg,this->pelec->f_en);
-        elecstate::ElecStateLCAO<std::complex<double>>* pelec_lcao
-            = dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec);
-        elecstate::cal_dm_psi(this->ParaV, pelec_lcao->wg, *psi_t, *(pelec_lcao->get_DM()));
-        pelec_lcao->get_DM()->cal_DMR();
+
+        elecstate::cal_dm_psi(this->ParaV, this->pelec->wg, *psi_t, *this->dm_);
+
+        this->dm_->cal_DMR();
+
         this->cal_mi_lcao(i_step);
     }
     else
@@ -431,27 +434,27 @@ void spinconstrain::SpinConstrain<std::complex<double>>::update_psi_charge(const
 
             if(pw_solve)
             {
-                hsolver::HSolverPW<std::complex<double>, base_device::DEVICE_CPU> hsolver_pw_obj(this->pw_wfc_,
-                                                 PARAM.inp.calculation,
-                                                 PARAM.inp.basis_type,
-                                                 PARAM.inp.ks_solver,
-                                                 false,
-                                                 PARAM.globalv.use_uspp,
-                                                 PARAM.inp.nspin,
-                                                 hsolver::DiagoIterAssist<std::complex<double>, base_device::DEVICE_CPU>::SCF_ITER,
-                                                 hsolver::DiagoIterAssist<std::complex<double>, base_device::DEVICE_CPU>::PW_DIAG_NMAX,
-                                                 hsolver::DiagoIterAssist<std::complex<double>, base_device::DEVICE_CPU>::PW_DIAG_THR,
-                                                 hsolver::DiagoIterAssist<std::complex<double>, base_device::DEVICE_CPU>::need_subspace);
+				hsolver::HSolverPW<std::complex<double>, base_device::DEVICE_CPU> hsolver_pw_obj(this->pw_wfc_,
+						PARAM.inp.calculation,
+						PARAM.inp.basis_type,
+						PARAM.inp.ks_solver,
+						false,
+						PARAM.globalv.use_uspp,
+						PARAM.inp.nspin,
+						hsolver::DiagoIterAssist<std::complex<double>, base_device::DEVICE_CPU>::SCF_ITER,
+						hsolver::DiagoIterAssist<std::complex<double>, base_device::DEVICE_CPU>::PW_DIAG_NMAX,
+						hsolver::DiagoIterAssist<std::complex<double>, base_device::DEVICE_CPU>::PW_DIAG_THR,
+						hsolver::DiagoIterAssist<std::complex<double>, base_device::DEVICE_CPU>::need_subspace);
 
-                hsolver_pw_obj.solve(hamilt_t,
-                         psi_t[0],
-                         this->pelec,
-                         this->pelec->ekb.c,
-                         GlobalV::RANK_IN_POOL,
-                         GlobalV::NPROC_IN_POOL,
-                         false,
-                         this->tpiba,
-                         this->get_nat());
+				hsolver_pw_obj.solve(hamilt_t,
+						psi_t[0],
+						this->pelec,
+						this->pelec->ekb.c,
+						GlobalV::RANK_IN_POOL,
+						GlobalV::NPROC_IN_POOL,
+						false,
+						this->tpiba,
+						this->get_nat());
             }
             else
             {// update charge density only
@@ -461,17 +464,17 @@ void spinconstrain::SpinConstrain<std::complex<double>>::update_psi_charge(const
 #if ((defined __CUDA) || (defined __ROCM))
         else
         {
-            base_device::DEVICE_GPU* ctx = {};
-            base_device::DEVICE_CPU* cpu_ctx = {};
-            psi::Psi<std::complex<double>, base_device::DEVICE_GPU>* psi_t = static_cast<psi::Psi<std::complex<double>, base_device::DEVICE_GPU>*>(this->psi);
-            hamilt::Hamilt<std::complex<double>, base_device::DEVICE_GPU>* hamilt_t = static_cast<hamilt::Hamilt<std::complex<double>, base_device::DEVICE_GPU>*>(this->p_hamilt);
-            auto* onsite_p = projectors::OnsiteProjector<double, base_device::DEVICE_GPU>::get_instance();
-            nbands = psi_t->get_nbands();
-            npol = psi_t->get_npol();
-            nkb = onsite_p->get_tot_nproj();
-            nk = psi_t->get_nk();
-            nh_iat = &onsite_p->get_nh(0);
-            size_becp = nbands * nkb * npol;
+			base_device::DEVICE_GPU* ctx = {};
+			base_device::DEVICE_CPU* cpu_ctx = {};
+			psi::Psi<std::complex<double>, base_device::DEVICE_GPU>* psi_t = static_cast<psi::Psi<std::complex<double>, base_device::DEVICE_GPU>*>(this->psi);
+			hamilt::Hamilt<std::complex<double>, base_device::DEVICE_GPU>* hamilt_t = static_cast<hamilt::Hamilt<std::complex<double>, base_device::DEVICE_GPU>*>(this->p_hamilt);
+			auto* onsite_p = projectors::OnsiteProjector<double, base_device::DEVICE_GPU>::get_instance();
+			nbands = psi_t->get_nbands();
+			npol = psi_t->get_npol();
+			nkb = onsite_p->get_tot_nproj();
+			nk = psi_t->get_nk();
+			nh_iat = &onsite_p->get_nh(0);
+			size_becp = nbands * nkb * npol;
 
             std::complex<double>* h_tmp = nullptr;
             std::complex<double>* s_tmp = nullptr;
diff --git a/source/source_lcao/module_deltaspin/init_sc.cpp b/source/source_lcao/module_deltaspin/init_sc.cpp
index ab192fe133..ac56047173 100644
--- a/source/source_lcao/module_deltaspin/init_sc.cpp
+++ b/source/source_lcao/module_deltaspin/init_sc.cpp
@@ -1,21 +1,24 @@
 #include "spin_constrain.h"
 
 // init sc
-template <typename FPTYPE>
-void spinconstrain::SpinConstrain<FPTYPE>::init_sc(double sc_thr_in,
-                                            int nsc_in,
-                                            int nsc_min_in,
-                                            double alpha_trial_in,
-                                            double sccut_in,
-                                            double sc_drop_thr_in,
-                                            const UnitCell& ucell,
-                                            Parallel_Orbitals* ParaV_in,
-                                            int nspin_in,
-                                            const K_Vectors& kv_in,
-                                            void* p_hamilt_in,
-                                            void* psi_in,
-                                            elecstate::ElecState* pelec_in,
-                                            ModulePW::PW_Basis_K* pw_wfc_in)
+template <typename TK>
+void spinconstrain::SpinConstrain<TK>::init_sc(double sc_thr_in,
+		int nsc_in,
+		int nsc_min_in,
+		double alpha_trial_in,
+		double sccut_in,
+		double sc_drop_thr_in,
+		const UnitCell& ucell,
+		Parallel_Orbitals* ParaV_in,
+		int nspin_in,
+		const K_Vectors& kv_in,
+		void* p_hamilt_in,
+		void* psi_in,
+#ifdef __LCAO
+		elecstate::DensityMatrix<TK, double>* dm_in, // mohan add 2025-11-03
+#endif
+		elecstate::ElecState* pelec_in,
+		ModulePW::PW_Basis_K* pw_wfc_in)
 {
     this->set_input_parameters(sc_thr_in, nsc_in, nsc_min_in, alpha_trial_in, sccut_in, sc_drop_thr_in);
     this->set_atomCounts(ucell.get_atom_Counts());
@@ -31,6 +34,9 @@ void spinconstrain::SpinConstrain<FPTYPE>::init_sc(double sc_thr_in,
     this->set_decay_grad();
     if(ParaV_in != nullptr) this->set_ParaV(ParaV_in);
     this->set_solver_parameters(kv_in, p_hamilt_in, psi_in, pelec_in);
+#ifdef __LCAO
+    this->dm_ = dm_in; // mohan add 2025-11-03
+#endif
 }
 
 template class spinconstrain::SpinConstrain<std::complex<double>>;
diff --git a/source/source_lcao/module_deltaspin/lambda_loop.cpp b/source/source_lcao/module_deltaspin/lambda_loop.cpp
index ac4fe88c9a..8b87f7eedd 100644
--- a/source/source_lcao/module_deltaspin/lambda_loop.cpp
+++ b/source/source_lcao/module_deltaspin/lambda_loop.cpp
@@ -99,7 +99,9 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
 
 
 template <>
-void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int outer_step, bool rerun)
+void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(
+        int outer_step,
+		bool rerun)
 {
     // init controlling parameters
     int nat = this->get_nat();
@@ -139,6 +141,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
         double duration = 0.0;
         if (i_step == -1)
         {
+
             this->cal_mw_from_lambda(i_step);
             spin = this->Mi_;
             where_fill_scalar_else_2d(this->constrain_, 0, zero, this->lambda_, initial_lambda);
@@ -151,7 +154,9 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
         {
             where_fill_scalar_else_2d(this->constrain_, 0, zero, delta_lambda, delta_lambda);
             add_scalar_multiply_2d(initial_lambda, delta_lambda, one, this->lambda_);
-            this->cal_mw_from_lambda(i_step, delta_lambda.data());
+
+            this->cal_mw_from_lambda(i_step);
+
             new_spin = this->Mi_;
             bool GradLessThanBound = this->check_gradient_decay(new_spin, spin, delta_lambda, dnu_last_step);
             if (i_step >= this->nsc_min_ && GradLessThanBound)
@@ -246,6 +251,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
 
         where_fill_scalar_else_2d(this->constrain_, 0, zero, delta_lambda, delta_lambda);
         add_scalar_multiply_2d(initial_lambda, delta_lambda, one, this->lambda_);
+
         this->cal_mw_from_lambda(i_step, delta_lambda.data());
 
         spin_plus = this->Mi_;
diff --git a/source/source_lcao/module_deltaspin/spin_constrain.cpp b/source/source_lcao/module_deltaspin/spin_constrain.cpp
index 6213f528f7..6b49487975 100644
--- a/source/source_lcao/module_deltaspin/spin_constrain.cpp
+++ b/source/source_lcao/module_deltaspin/spin_constrain.cpp
@@ -8,15 +8,15 @@
 namespace spinconstrain
 {
 
-template <typename FPTYPE>
-SpinConstrain<FPTYPE>& SpinConstrain<FPTYPE>::getScInstance()
+template <typename TK>
+SpinConstrain<TK>& SpinConstrain<TK>::getScInstance()
 {
-    static SpinConstrain<FPTYPE> instance; // Guaranteed to be created and destroyed only once
+    static SpinConstrain<TK> instance; // Guaranteed to be created and destroyed only once
     return instance;
 }
 
-template <typename FPTYPE>
-double SpinConstrain<FPTYPE>::cal_escon()
+template <typename TK>
+double SpinConstrain<TK>::cal_escon()
 {
     this->escon_ = 0.0;
     if (!this->is_Mi_converged)
@@ -33,30 +33,30 @@ double SpinConstrain<FPTYPE>::cal_escon()
     return this->escon_;
 }
 
-template <typename FPTYPE>
-double SpinConstrain<FPTYPE>::get_escon()
+template <typename TK>
+double SpinConstrain<TK>::get_escon()
 {
     return this->escon_;
 }
 
 // set atomCounts
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_atomCounts(const std::map<int, int>& atomCounts_in)
+template <typename TK>
+void SpinConstrain<TK>::set_atomCounts(const std::map<int, int>& atomCounts_in)
 {
     this->atomCounts.clear();
     this->atomCounts = atomCounts_in;
 }
 
 // get atomCounts
-template <typename FPTYPE>
-const std::map<int, int>& SpinConstrain<FPTYPE>::get_atomCounts() const
+template <typename TK>
+const std::map<int, int>& SpinConstrain<TK>::get_atomCounts() const
 {
     return this->atomCounts;
 }
 
 /// set nspin
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_nspin(int nspin_in)
+template <typename TK>
+void SpinConstrain<TK>::set_nspin(int nspin_in)
 {
     if (nspin_in != 4 && nspin_in != 2)
     {
@@ -66,14 +66,14 @@ void SpinConstrain<FPTYPE>::set_nspin(int nspin_in)
 }
 
 /// get nspin
-template <typename FPTYPE>
-int SpinConstrain<FPTYPE>::get_nspin()
+template <typename TK>
+int SpinConstrain<TK>::get_nspin()
 {
     return this->nspin_;
 }
 
-template <typename FPTYPE>
-int SpinConstrain<FPTYPE>::get_nat()
+template <typename TK>
+int SpinConstrain<TK>::get_nat()
 {
     int nat = 0;
     for (std::map<int, int>::iterator it = this->atomCounts.begin(); it != this->atomCounts.end(); ++it)
@@ -83,14 +83,14 @@ int SpinConstrain<FPTYPE>::get_nat()
     return nat;
 }
 
-template <typename FPTYPE>
-int SpinConstrain<FPTYPE>::get_ntype()
+template <typename TK>
+int SpinConstrain<TK>::get_ntype()
 {
     return this->atomCounts.size();
 }
 
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::check_atomCounts()
+template <typename TK>
+void SpinConstrain<TK>::check_atomCounts()
 {
     if (!this->atomCounts.size())
     {
@@ -116,8 +116,8 @@ void SpinConstrain<FPTYPE>::check_atomCounts()
 }
 
 // get iat
-template <typename FPTYPE>
-int SpinConstrain<FPTYPE>::get_iat(int itype, int atom_index)
+template <typename TK>
+int SpinConstrain<TK>::get_iat(int itype, int atom_index)
 {
     if (itype < 0 || itype >= this->get_ntype())
     {
@@ -141,38 +141,38 @@ int SpinConstrain<FPTYPE>::get_iat(int itype, int atom_index)
 }
 
 // set orbitalCounts
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_orbitalCounts(const std::map<int, int>& orbitalCounts_in)
+template <typename TK>
+void SpinConstrain<TK>::set_orbitalCounts(const std::map<int, int>& orbitalCounts_in)
 {
     this->orbitalCounts.clear();
     this->orbitalCounts = orbitalCounts_in;
 }
 
 // get orbitalCounts
-template <typename FPTYPE>
-const std::map<int, int>& SpinConstrain<FPTYPE>::get_orbitalCounts() const
+template <typename TK>
+const std::map<int, int>& SpinConstrain<TK>::get_orbitalCounts() const
 {
     return this->orbitalCounts;
 }
 
 // set lnchiCounts
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_lnchiCounts(const std::map<int, std::map<int, int>>& lnchiCounts_in)
+template <typename TK>
+void SpinConstrain<TK>::set_lnchiCounts(const std::map<int, std::map<int, int>>& lnchiCounts_in)
 {
     this->lnchiCounts.clear();
     this->lnchiCounts = lnchiCounts_in;
 }
 
 // get lnchiCounts
-template <typename FPTYPE>
-const std::map<int, std::map<int, int>>& SpinConstrain<FPTYPE>::get_lnchiCounts() const
+template <typename TK>
+const std::map<int, std::map<int, int>>& SpinConstrain<TK>::get_lnchiCounts() const
 {
     return this->lnchiCounts;
 }
 
 // set sc_lambda from ScData
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_sc_lambda()
+template <typename TK>
+void SpinConstrain<TK>::set_sc_lambda()
 {
     this->check_atomCounts();
     int nat = this->get_nat();
@@ -194,8 +194,8 @@ void SpinConstrain<FPTYPE>::set_sc_lambda()
 }
 
 // set target_mag from ScData
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_target_mag()
+template <typename TK>
+void SpinConstrain<TK>::set_target_mag()
 {
     this->check_atomCounts();
     int nat = this->get_nat();
@@ -234,8 +234,8 @@ void SpinConstrain<FPTYPE>::set_target_mag()
 }
 
 // set constrain from ScData
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_constrain()
+template <typename TK>
+void SpinConstrain<TK>::set_constrain()
 {
     this->check_atomCounts();
     int nat = this->get_nat();
@@ -266,8 +266,8 @@ void SpinConstrain<FPTYPE>::set_constrain()
 }
 
 // set sc_lambda from variable
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_sc_lambda(const ModuleBase::Vector3<double>* lambda_in, int nat_in)
+template <typename TK>
+void SpinConstrain<TK>::set_sc_lambda(const ModuleBase::Vector3<double>* lambda_in, int nat_in)
 {
     this->check_atomCounts();
     int nat = this->get_nat();
@@ -283,8 +283,8 @@ void SpinConstrain<FPTYPE>::set_sc_lambda(const ModuleBase::Vector3<double>* lam
 }
 
 // set target_mag from variable
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_target_mag(const ModuleBase::Vector3<double>* target_mag_in, int nat_in)
+template <typename TK>
+void SpinConstrain<TK>::set_target_mag(const ModuleBase::Vector3<double>* target_mag_in, int nat_in)
 {
     this->check_atomCounts();
     int nat = this->get_nat();
@@ -299,8 +299,8 @@ void SpinConstrain<FPTYPE>::set_target_mag(const ModuleBase::Vector3<double>* ta
     }
 }
 
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_target_mag(const std::vector<ModuleBase::Vector3<double>>& target_mag_in)
+template <typename TK>
+void SpinConstrain<TK>::set_target_mag(const std::vector<ModuleBase::Vector3<double>>& target_mag_in)
 {
     int nat = this->get_nat();
     assert(target_mag_in.size() == nat);
@@ -324,8 +324,8 @@ void SpinConstrain<FPTYPE>::set_target_mag(const std::vector<ModuleBase::Vector3
 }
 
 /// set constrain from variable
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_constrain(const ModuleBase::Vector3<int>* constrain_in, int nat_in)
+template <typename TK>
+void SpinConstrain<TK>::set_constrain(const ModuleBase::Vector3<int>* constrain_in, int nat_in)
 {
     this->check_atomCounts();
     int nat = this->get_nat();
@@ -340,28 +340,28 @@ void SpinConstrain<FPTYPE>::set_constrain(const ModuleBase::Vector3<int>* constr
     }
 }
 
-template <typename FPTYPE>
-const std::vector<ModuleBase::Vector3<double>>& SpinConstrain<FPTYPE>::get_sc_lambda() const
+template <typename TK>
+const std::vector<ModuleBase::Vector3<double>>& SpinConstrain<TK>::get_sc_lambda() const
 {
     return this->lambda_;
 }
 
-template <typename FPTYPE>
-const std::vector<ModuleBase::Vector3<double>>& SpinConstrain<FPTYPE>::get_target_mag() const
+template <typename TK>
+const std::vector<ModuleBase::Vector3<double>>& SpinConstrain<TK>::get_target_mag() const
 {
     return this->target_mag_;
 }
 
 /// get_constrain
-template <typename FPTYPE>
-const std::vector<ModuleBase::Vector3<int>>& SpinConstrain<FPTYPE>::get_constrain() const
+template <typename TK>
+const std::vector<ModuleBase::Vector3<int>>& SpinConstrain<TK>::get_constrain() const
 {
     return this->constrain_;
 }
 
 /// zero atomic magnetic moment
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::zero_Mi()
+template <typename TK>
+void SpinConstrain<TK>::zero_Mi()
 {
     this->check_atomCounts();
     int nat = this->get_nat();
@@ -377,15 +377,15 @@ void SpinConstrain<FPTYPE>::zero_Mi()
 /// get grad_decay
 /// this function can only be called by the root process because only
 /// root process reads the ScDecayGrad from json file
-template <typename FPTYPE>
-double SpinConstrain<FPTYPE>::get_decay_grad(int itype)
+template <typename TK>
+double SpinConstrain<TK>::get_decay_grad(int itype)
 {
     return this->ScDecayGrad[itype];
 }
 
 /// set grad_decy
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_decay_grad()
+template <typename TK>
+void SpinConstrain<TK>::set_decay_grad()
 {
     this->check_atomCounts();
     int ntype = this->get_ntype();
@@ -397,15 +397,15 @@ void SpinConstrain<FPTYPE>::set_decay_grad()
 }
 
 /// get decay_grad
-template <typename FPTYPE>
-const std::vector<double>& SpinConstrain<FPTYPE>::get_decay_grad()
+template <typename TK>
+const std::vector<double>& SpinConstrain<TK>::get_decay_grad()
 {
     return this->decay_grad_;
 }
 
 /// set grad_decy from variable
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_decay_grad(const double* decay_grad_in, int ntype_in)
+template <typename TK>
+void SpinConstrain<TK>::set_decay_grad(const double* decay_grad_in, int ntype_in)
 {
     this->check_atomCounts();
     int ntype = this->get_ntype();
@@ -421,8 +421,8 @@ void SpinConstrain<FPTYPE>::set_decay_grad(const double* decay_grad_in, int ntyp
 }
 
 /// @brief  set input parameters
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_input_parameters(double sc_thr_in,
+template <typename TK>
+void SpinConstrain<TK>::set_input_parameters(double sc_thr_in,
                                                  int nsc_in,
                                                  int nsc_min_in,
                                                  double alpha_trial_in,
@@ -438,56 +438,56 @@ void SpinConstrain<FPTYPE>::set_input_parameters(double sc_thr_in,
 }
 
 /// get sc_thr
-template <typename FPTYPE>
-double SpinConstrain<FPTYPE>::get_sc_thr()
+template <typename TK>
+double SpinConstrain<TK>::get_sc_thr()
 {
     return this->sc_thr_;
 }
 
 /// get nsc
-template <typename FPTYPE>
-int SpinConstrain<FPTYPE>::get_nsc()
+template <typename TK>
+int SpinConstrain<TK>::get_nsc()
 {
     return this->nsc_;
 }
 
 /// get nsc_min
-template <typename FPTYPE>
-int SpinConstrain<FPTYPE>::get_nsc_min()
+template <typename TK>
+int SpinConstrain<TK>::get_nsc_min()
 {
     return this->nsc_min_;
 }
 
 /// get alpha_trial
-template <typename FPTYPE>
-double SpinConstrain<FPTYPE>::get_alpha_trial()
+template <typename TK>
+double SpinConstrain<TK>::get_alpha_trial()
 {
     return this->alpha_trial_;
 }
 
 /// get sccut
-template <typename FPTYPE>
-double SpinConstrain<FPTYPE>::get_sccut()
+template <typename TK>
+double SpinConstrain<TK>::get_sccut()
 {
     return this->restrict_current_;
 }
 
 /// set sc_drop_thr
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_sc_drop_thr(double sc_drop_thr_in)
+template <typename TK>
+void SpinConstrain<TK>::set_sc_drop_thr(double sc_drop_thr_in)
 {
     this->sc_drop_thr_ = sc_drop_thr_in;
 }
 
 /// get sc_drop_thr
-template <typename FPTYPE>
-double SpinConstrain<FPTYPE>::get_sc_drop_thr()
+template <typename TK>
+double SpinConstrain<TK>::get_sc_drop_thr()
 {
     return this->sc_drop_thr_;
 }
 
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_solver_parameters(const K_Vectors& kv_in,
+template <typename TK>
+void SpinConstrain<TK>::set_solver_parameters(const K_Vectors& kv_in,
                                                   void* p_hamilt_in,
                                                   void* psi_in,
                                                   elecstate::ElecState* pelec_in)
@@ -499,8 +499,8 @@ void SpinConstrain<FPTYPE>::set_solver_parameters(const K_Vectors& kv_in,
 }
 
 /// @brief  set ParaV
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::set_ParaV(Parallel_Orbitals* ParaV_in)
+template <typename TK>
+void SpinConstrain<TK>::set_ParaV(Parallel_Orbitals* ParaV_in)
 {
     this->ParaV = ParaV_in;
     int nloc = this->ParaV->nloc;
@@ -511,8 +511,8 @@ void SpinConstrain<FPTYPE>::set_ParaV(Parallel_Orbitals* ParaV_in)
 }
 
 /// print Mi
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::print_Mi(std::ofstream& ofs_running)
+template <typename TK>
+void SpinConstrain<TK>::print_Mi(std::ofstream& ofs_running)
 {
     this->check_atomCounts();
     int nat = this->get_nat();
@@ -556,8 +556,8 @@ void SpinConstrain<FPTYPE>::print_Mi(std::ofstream& ofs_running)
 }
 
 /// print magnetic force (defined as \frac{\delta{L}}/{\delta{Mi}} = -lambda[iat])
-template <typename FPTYPE>
-void SpinConstrain<FPTYPE>::print_Mag_Force(std::ofstream& ofs_running)
+template <typename TK>
+void SpinConstrain<TK>::print_Mag_Force(std::ofstream& ofs_running)
 {
     this->check_atomCounts();
     int nat = this->get_nat();
diff --git a/source/source_lcao/module_deltaspin/spin_constrain.h b/source/source_lcao/module_deltaspin/spin_constrain.h
index bcc4ed9ec1..f4868aa93a 100644
--- a/source/source_lcao/module_deltaspin/spin_constrain.h
+++ b/source/source_lcao/module_deltaspin/spin_constrain.h
@@ -14,12 +14,16 @@
 #include "source_hamilt/operator.h"
 #include "source_estate/elecstate.h"
 
+#ifdef __LCAO
+#include "source_estate/module_dm/density_matrix.h" // mohan add 2025-11-02
+#endif
+
 namespace spinconstrain
 {
 
 struct ScAtomData;
 
-template <typename FPTYPE>
+template <typename TK>
 class SpinConstrain
 {
 public:
@@ -39,7 +43,10 @@ class SpinConstrain
                const K_Vectors& kv_in,
                void* p_hamilt_in,
                void* psi_in,
-               elecstate::ElecState* pelec_in,
+#ifdef __LCAO
+			   elecstate::DensityMatrix<TK, double> *dm_in, // mohan add 2025-11-02
+#endif
+			   elecstate::ElecState* pelec_in,
                ModulePW::PW_Basis_K* pw_wfc_in = nullptr);
 
   /// @brief calculate the magnetization of each atom with real space projection method for LCAO base
@@ -49,7 +56,8 @@ class SpinConstrain
 
   void cal_mi_pw();
 
-  void cal_mw_from_lambda(int i_step, const ModuleBase::Vector3<double>* delta_lambda = nullptr);
+  void cal_mw_from_lambda(int i_step, 
+		  const ModuleBase::Vector3<double>* delta_lambda = nullptr);
 
   /**
    * @brief calculate the energy of \sum_i \lambda_i * Mi
@@ -60,13 +68,17 @@ class SpinConstrain
 
   double get_escon();
 
-  void run_lambda_loop(int outer_step, bool rerun = true);
+  void run_lambda_loop(int outer_step, 
+		  bool rerun = true);
 
   /// @brief update the charge density for LCAO base with new lambda
   /// update the charge density and psi for PW base with new lambda
   void update_psi_charge(const ModuleBase::Vector3<double>* delta_lambda, bool pw_solve = true);
 
-  void calculate_delta_hcc(std::complex<double>* h_tmp, const std::complex<double>* becp_k, const ModuleBase::Vector3<double>* delta_lambda, const int nbands, const int nkb, const int* nh_iat);
+  void calculate_delta_hcc(std::complex<double>* h_tmp, 
+		  const std::complex<double>* becp_k, 
+		  const ModuleBase::Vector3<double>* delta_lambda, 
+		  const int nbands, const int nkb, const int* nh_iat);
 
   /// lambda loop helper functions
   bool check_rms_stop(int outer_step, int i_step, double rms_error, double duration, double total_duration);
@@ -109,6 +121,9 @@ class SpinConstrain
     void* psi = nullptr;
     elecstate::ElecState* pelec = nullptr;
     ModulePW::PW_Basis_K* pw_wfc_ = nullptr;
+#ifdef __LCAO
+    elecstate::DensityMatrix<TK, double>* dm_;
+#endif
     double tpiba = 0.0; /// save ucell.tpiba
     const double meV_to_Ry = 7.349864435130999e-05;
     K_Vectors kv_;
@@ -240,20 +255,20 @@ class SpinConstrain
   public:
     /// @brief save operator for spin-constrained DFT
     /// @param op_in the base pointer of operator, actual type should be DeltaSpin<OperatorLCAO<TK, TR>>*
-    void set_operator(hamilt::Operator<FPTYPE>* op_in);
+    void set_operator(hamilt::Operator<TK>* op_in);
     /// @brief set is_Mi_converged
     void set_mag_converged(bool is_Mi_converged_in){this->is_Mi_converged = is_Mi_converged_in;}
     /// @brief get is_Mi_converged
     bool mag_converged() const {return this->is_Mi_converged;}
   private:
     /// operator for spin-constrained DFT, used for calculating current atomic magnetic moment
-    hamilt::Operator<FPTYPE>* p_operator = nullptr;
+    hamilt::Operator<TK>* p_operator = nullptr;
     /// @brief if atomic magnetic moment is converged
     bool is_Mi_converged = false;
 
-    FPTYPE* sub_h_save;
-    FPTYPE* sub_s_save;
-    FPTYPE* becp_save;
+    TK* sub_h_save;
+    TK* sub_s_save;
+    TK* becp_save;
 };
 
 
diff --git a/source/source_lcao/module_deltaspin/template_helpers.cpp b/source/source_lcao/module_deltaspin/template_helpers.cpp
index 437ecec662..83e5f17f75 100644
--- a/source/source_lcao/module_deltaspin/template_helpers.cpp
+++ b/source/source_lcao/module_deltaspin/template_helpers.cpp
@@ -1,7 +1,8 @@
 #include "spin_constrain.h"
 
 template <>
-void spinconstrain::SpinConstrain<double>::cal_mw_from_lambda(int i_step, const ModuleBase::Vector3<double>* delta_lambda)
+void spinconstrain::SpinConstrain<double>::cal_mw_from_lambda(int i_step, 
+		const ModuleBase::Vector3<double>* delta_lambda)
 {
 }
 
@@ -11,7 +12,8 @@ void spinconstrain::SpinConstrain<double>::cal_mi_lcao(const int& step, bool pri
 }
 
 template <>
-void spinconstrain::SpinConstrain<double>::run_lambda_loop(int outer_step, bool rerun)
+void spinconstrain::SpinConstrain<double>::run_lambda_loop(int outer_step, 
+		bool rerun)
 {
 }
 
@@ -60,4 +62,4 @@ bool spinconstrain::SpinConstrain<double>::check_gradient_decay(
     bool print)
 {
     return false;
-}
\ No newline at end of file
+}
diff --git a/source/source_lcao/module_dftu/dftu.h b/source/source_lcao/module_dftu/dftu.h
index b269acc665..457debc8bd 100644
--- a/source/source_lcao/module_dftu/dftu.h
+++ b/source/source_lcao/module_dftu/dftu.h
@@ -7,7 +7,7 @@
 #ifdef __LCAO
 #include "source_estate/module_charge/charge_mixing.h"
 #include "source_hamilt/hamilt.h"
-#include "source_estate/elecstate.h"
+// #include "source_estate/elecstate.h" // mohan update 2025-11-02
 #include "source_lcao/module_hcontainer/hcontainer.h"
 #include "source_estate/module_dm/density_matrix.h"
 #include "source_lcao/force_stress_arrays.h" // mohan add 2024-06-15
@@ -200,8 +200,9 @@ class DFTU
  public:
    void force_stress(const UnitCell& ucell,
                      const Grid_Driver& gd,
-                     const elecstate::ElecState* pelec,
-                     const Parallel_Orbitals& pv,
+					 std::vector<std::vector<double>>* dmk_d, // mohan modify 2025-11-02
+					 std::vector<std::vector<std::complex<double>>>* dmk_c, // dmat.get_dm()->get_DMK_vector();
+					 const Parallel_Orbitals& pv,
                      ForceStressArrays& fsr,
                      ModuleBase::matrix& force_dftu,
                      ModuleBase::matrix& stress_dftu,
diff --git a/source/source_lcao/module_dftu/dftu_force.cpp b/source/source_lcao/module_dftu/dftu_force.cpp
index 5f08b8b545..5d76611470 100644
--- a/source/source_lcao/module_dftu/dftu_force.cpp
+++ b/source/source_lcao/module_dftu/dftu_force.cpp
@@ -1,8 +1,5 @@
-//==========================================================
-// Author:Xin Qu
 #include "source_io/module_parameter/parameter.h"
-// DATE : 2019-12-10
-//==========================================================
+
 #ifdef __LCAO
 #include "dftu.h"
 #include "source_base/constants.h"
@@ -24,6 +21,7 @@
 #include <stdio.h>
 #include <string.h>
 
+
 extern "C"
 {
     // I'm not sure what's happenig here, but the interface in scalapack_connecter.h
@@ -74,7 +72,8 @@ namespace ModuleDFTU
 
 void DFTU::force_stress(const UnitCell& ucell,
                         const Grid_Driver& gd,
-                        const elecstate::ElecState* pelec,
+                        std::vector<std::vector<double>>* dmk_d, // mohan modify 2025-11-02
+                        std::vector<std::vector<std::complex<double>>>* dmk_c, // dmat.get_dm()->get_DMK_vector();
                         const Parallel_Orbitals& pv,
                         ForceStressArrays& fsr, // mohan add 2024-06-16
                         ModuleBase::matrix& force_dftu,
@@ -114,29 +113,12 @@ void DFTU::force_stress(const UnitCell& ucell,
 
             this->cal_VU_pot_mat_real(spin, false, VU);
 
-            const std::vector<std::vector<double>>& dmk
-                = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec)->get_DM()->get_DMK_vector();
-
 #ifdef __MPI
-            pdgemm_(&transT,
-                    &transN,
-                    &nlocal,
-                    &nlocal,
-                    &nlocal,
-                    &alpha,
-                    dmk[spin].data(),
-                    &one_int,
-                    &one_int,
-                    pv.desc,
-                    VU,
-                    &one_int,
-                    &one_int,
-                    pv.desc,
-                    &beta,
-                    &rho_VU[0],
-                    &one_int,
-                    &one_int,
-                    pv.desc);
+            pdgemm_(&transT, &transN, &nlocal, &nlocal, &nlocal,
+                    &alpha, (*dmk_d)[spin].data(), &one_int, &one_int, // important to add () outside *dmk_d, mohan note 20251103
+                    pv.desc, VU, &one_int, &one_int,
+                    pv.desc, &beta, &rho_VU[0],
+                    &one_int, &one_int, pv.desc);
 #endif
 
             delete[] VU;
@@ -178,31 +160,12 @@ void DFTU::force_stress(const UnitCell& ucell,
 
             this->cal_VU_pot_mat_complex(spin, false, VU);
 
-            const std::vector<std::vector<std::complex<double>>>& dmk
-                = dynamic_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(pelec)
-                      ->get_DM()
-                      ->get_DMK_vector();
 
 #ifdef __MPI
-            pzgemm_(&transT,
-                    &transN,
-                    &nlocal,
-                    &nlocal,
-                    &nlocal,
-                    &alpha,
-                    dmk[ik].data(),
-                    &one_int,
-                    &one_int,
-                    pv.desc,
-                    VU,
-                    &one_int,
-                    &one_int,
-                    pv.desc,
-                    &beta,
-                    &rho_VU[0],
-                    &one_int,
-                    &one_int,
-                    pv.desc);
+            pzgemm_(&transT, &transN, &nlocal, &nlocal, &nlocal,
+                    &alpha, (*dmk_c)[ik].data(), &one_int, &one_int, // important to add (), 20251103
+                    pv.desc, VU, &one_int, &one_int, pv.desc, &beta,
+                    &rho_VU[0], &one_int, &one_int, pv.desc);
 #endif
 
             delete[] VU;
diff --git a/source/source_lcao/module_ri/Exx_LRI_interface.h b/source/source_lcao/module_ri/Exx_LRI_interface.h
index ff76dbbc7d..15302677b5 100644
--- a/source/source_lcao/module_ri/Exx_LRI_interface.h
+++ b/source/source_lcao/module_ri/Exx_LRI_interface.h
@@ -4,6 +4,7 @@
 #include "Exx_LRI.h"
 #include "source_lcao/module_ri/Mix_DMk_2D.h"
 #include "source_lcao/module_ri/module_exx_symmetry/symmetry_rotation.h"
+#include "source_estate/module_dm/density_matrix.h" // mohan add 2025-11-04
 #include <memory>
 
 class LCAO_Matrix;
@@ -96,7 +97,8 @@ class Exx_LRI_Interface
     void exx_iter_finish(const K_Vectors& kv,
                          const UnitCell& ucell,
                          hamilt::Hamilt<T>& hamilt,
-                         elecstate::ElecState& elec,
+						 elecstate::ElecState& elec,
+						 elecstate::DensityMatrix<T,double>* dm, // mohan add 2025-11-04
                          Charge_Mixing& chgmix,
                          const double& scf_ene_thr,
                          int& iter,
diff --git a/source/source_lcao/module_ri/Exx_LRI_interface.hpp b/source/source_lcao/module_ri/Exx_LRI_interface.hpp
index 89db3e5417..802f02a4f8 100644
--- a/source/source_lcao/module_ri/Exx_LRI_interface.hpp
+++ b/source/source_lcao/module_ri/Exx_LRI_interface.hpp
@@ -262,14 +262,15 @@ void Exx_LRI_Interface<T, Tdata>::exx_hamilt2rho(elecstate::ElecState& elec, con
 
 template<typename T, typename Tdata>
 void Exx_LRI_Interface<T, Tdata>::exx_iter_finish(const K_Vectors& kv,
-                                                  const UnitCell& ucell,
-                                                  hamilt::Hamilt<T>& hamilt,
-                                                  elecstate::ElecState& elec,
-                                                  Charge_Mixing& chgmix,
-                                                  const double& scf_ene_thr,
-                                                  int& iter,
-                                                  const int istep,
-                                                  bool& conv_esolver)
+		const UnitCell& ucell,
+		hamilt::Hamilt<T>& hamilt,
+		elecstate::ElecState& elec,
+		elecstate::DensityMatrix<T,double>* dm, // mohan add 2025-11-04
+		Charge_Mixing& chgmix,
+		const double& scf_ene_thr,
+		int& iter,
+		const int istep,
+		bool& conv_esolver)
 {
     ModuleBase::TITLE("Exx_LRI_Interface","exx_iter_finish");
     if (GlobalC::restart.info_save.save_H && (this->two_level_step > 0 || istep > 0)
@@ -312,11 +313,12 @@ void Exx_LRI_Interface<T, Tdata>::exx_iter_finish(const K_Vectors& kv,
         {
             chgmix.close_kerker_gg0();
         }
-        this->dm_last_step = dynamic_cast<const elecstate::ElecStateLCAO<T>*>(&elec)->get_DM();
+// mohan update 2025-11-04
+        this->dm_last_step = dm;
         conv_esolver = this->exx_after_converge(
             ucell,
             hamilt,
-            *dynamic_cast<const elecstate::ElecStateLCAO<T>*>(&elec)->get_DM(),
+            *dm,
             kv,
             PARAM.inp.nspin,
             iter,
diff --git a/source/source_lcao/setup_dm.cpp b/source/source_lcao/setup_dm.cpp
index 231c8b0615..a91b43ab96 100644
--- a/source/source_lcao/setup_dm.cpp
+++ b/source/source_lcao/setup_dm.cpp
@@ -1,6 +1,6 @@
 #include "source_lcao/setup_dm.h"
 
-#include "cal_dm.h"
+#include "source_estate/cal_dm.h"
 #include "source_base/timer.h"
 #include "source_estate/module_dm/cal_dm_psi.h"
 #include "source_hamilt/module_xc/xc_functional.h"
@@ -15,11 +15,12 @@
 namespace LCAO_domain
 {
 
+// change init_dm to allocate_dm, mohan 2025-10-31
 template <typename TK>
-void Setup_DM<TK>::init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin)
+void Setup_DM<TK>::allocate_dm(const K_Vectors* kv, const Parallel_Orbitals* pv, const int nspin)
 {
     const int nspin_dm = nspin == 2 ? 2 : 1;
-    this->dm = new DensityMatrix<TK, double>(paraV, nspin_dm, kv->kvec_d, kv->get_nks() / nspin_dm);
+    this->dm = new elecstate::DensityMatrix<TK, double>(pv, nspin_dm, kv->kvec_d, kv->get_nks() / nspin_dm);
 }
 
 template class Setup_DM<double>;               // Gamma_only case
diff --git a/source/source_lcao/setup_dm.h b/source/source_lcao/setup_dm.h
index 0e4da43497..672a50c878 100644
--- a/source/source_lcao/setup_dm.h
+++ b/source/source_lcao/setup_dm.h
@@ -1,11 +1,13 @@
 #ifndef SETUP_DM_H
 #define SETUP_DM_H 
 
+#include "source_cell/klist.h" 
+#include "source_basis/module_ao/parallel_orbitals.h"
 #include "source_estate/module_dm/density_matrix.h"
 
 #include <vector>
 
-namespace LCAO_DOMAIN
+namespace LCAO_domain
 {
 template <typename TK>
 class Setup_DM
@@ -18,23 +20,16 @@ class Setup_DM
 
     ~Setup_DM()
     {
-        if (this->DM != nullptr)
+        if (this->dm != nullptr)
         {
             delete this->dm;
         }
     }
 
-    // initial density matrix
-    void init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin);
+    // allocate density matrix
+    void allocate_dm(const K_Vectors* kv, const Parallel_Orbitals* pv, const int nspin);
 
-    DensityMatrix<TK, double>* get_dm() const
-    {
-        return const_cast<DensityMatrix<TK, double>*>(this->dm);
-    }
-
-    private:
-
-    DensityMatrix<TK, double>* dm = nullptr;
+    elecstate::DensityMatrix<TK, double>* dm = nullptr;
 
 };
 
diff --git a/source/source_psi/test/psi_initializer_unit_test.cpp b/source/source_psi/test/psi_initializer_unit_test.cpp
index 44d3a6eb9c..0bfad3d371 100644
--- a/source/source_psi/test/psi_initializer_unit_test.cpp
+++ b/source/source_psi/test/psi_initializer_unit_test.cpp
@@ -80,7 +80,7 @@ InfoNonlocal::~InfoNonlocal() {}
 #endif
 Structure_Factor::Structure_Factor() {}
 Structure_Factor::~Structure_Factor() {}
-void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const Parallel_Grid&, const ModulePW::PW_Basis* rho_basis) {}
+void Structure_Factor::setup(const UnitCell* Ucell, const Parallel_Grid&, const ModulePW::PW_Basis* rho_basis) {}
 std::complex<double>* Structure_Factor::get_sk(int ik, int it, int ia, ModulePW::PW_Basis_K const*wfc_basis) const
 {
     int npw = wfc_basis->npwk[ik];
@@ -546,4 +546,4 @@ int main(int argc, char** argv)
 #endif
 
     return result;
-}
\ No newline at end of file
+}
diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp
index 1eef712283..fb20b547e8 100644
--- a/source/source_pw/module_pwdft/setup_pot.cpp
+++ b/source/source_pw/module_pwdft/setup_pot.cpp
@@ -102,11 +102,14 @@ void pw::setup_pot(const int istep,
                    PARAM.inp.sccut,
                    PARAM.inp.sc_drop_thr,
                    ucell,
-                   nullptr,
+                   nullptr, // parallel orbitals
                    PARAM.inp.nspin,
                    kv,
                    p_hamilt,
                    kspw_psi,
+#ifdef __LCAO
+                   nullptr, // density matrix, not useful in LCAO, mohan note 2025-11-03
+#endif
                    pelec,
                    pw_wfc);
     }
diff --git a/source/source_pw/module_pwdft/structure_factor.cpp b/source/source_pw/module_pwdft/structure_factor.cpp
index b97bdee944..2823d11e36 100644
--- a/source/source_pw/module_pwdft/structure_factor.cpp
+++ b/source/source_pw/module_pwdft/structure_factor.cpp
@@ -46,7 +46,7 @@ Structure_Factor::~Structure_Factor()
 // called in input.cpp
 void Structure_Factor::set(const ModulePW::PW_Basis* rho_basis_in, const int& nbspline_in)
 {
-    ModuleBase::TITLE("PW_Basis","set");
+    ModuleBase::TITLE("Structure_Factor","set");
     this->rho_basis = rho_basis_in;
     this->nbspline = nbspline_in;
     return;
@@ -54,10 +54,11 @@ void Structure_Factor::set(const ModulePW::PW_Basis* rho_basis_in, const int& nb
 
 // Peize Lin optimize and add OpenMP 2021.04.01
 //  Calculate structure factor
-void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const Parallel_Grid& pgrid, const ModulePW::PW_Basis* rho_basis)
+void Structure_Factor::setup(const UnitCell* Ucell, const Parallel_Grid& pgrid, const ModulePW::PW_Basis* rho_basis)
 {
-    ModuleBase::TITLE("PW_Basis","setup_structure_factor");
-    ModuleBase::timer::tick("PW_Basis","setup_struc_factor");
+    ModuleBase::TITLE("Structure_Factor","setup");
+    ModuleBase::timer::tick("Structure_Factor","setup");
+
     const std::complex<double> ci_tpi = ModuleBase::NEG_IMAG_UNIT * ModuleBase::TWO_PI;
     this->ucell = Ucell;
     this->strucFac.create(Ucell->ntype, rho_basis->npw);
@@ -66,9 +67,15 @@ void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const Paral
 //	std::string outstr;
 //	outstr = PARAM.globalv.global_out_dir + "strucFac.dat"; 
 //	std::ofstream ofs( outstr.c_str() ) ;
-    bool usebspline;
-    if(nbspline > 0) {   usebspline = true;
-    } else {    usebspline = false;}
+	bool usebspline;
+	if(nbspline > 0) 
+	{   
+		usebspline = true;
+	} 
+	else 
+	{    
+		usebspline = false;
+	}
     
     if(usebspline)
     {
@@ -100,12 +107,15 @@ void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const Paral
 
 //	ofs.close();
 
-    int i,j; //ng;
+    int i=0;
+    int j=0;
+ 
     this->eigts1.create(Ucell->nat, 2*rho_basis->nx + 1);
     this->eigts2.create(Ucell->nat, 2*rho_basis->ny + 1);
     this->eigts3.create(Ucell->nat, 2*rho_basis->nz + 1);
 
-    ModuleBase::Memory::record("SF::eigts123",sizeof(std::complex<double>) * (Ucell->nat*2 * (rho_basis->nx + rho_basis->ny + rho_basis->nz) + 3));
+    ModuleBase::Memory::record("SF::eigts123",sizeof(std::complex<double>) 
+    * (Ucell->nat*2 * (rho_basis->nx + rho_basis->ny + rho_basis->nz) + 3));
 
     ModuleBase::Vector3<double> gtau;
     int inat = 0;
@@ -177,7 +187,7 @@ void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const Paral
         this->z_eigts3 = this->eigts3.c;
         // There's no need to delete double precision pointers while in a CPU environment.
     }
-    ModuleBase::timer::tick("PW_Basis","setup_struc_factor"); 
+    ModuleBase::timer::tick("Structure_Factor","setup");
     return;
 }
 
@@ -298,7 +308,7 @@ void Structure_Factor::bspline_sf(const int norder,
     return;
 }
 
-void Structure_Factor:: bsplinecoef(std::complex<double> *b1, std::complex<double> *b2, std::complex<double> *b3, 
+void Structure_Factor::bsplinecoef(std::complex<double> *b1, std::complex<double> *b2, std::complex<double> *b3, 
                         const int nx, const int ny, const int nz, const int norder)
 {
     const std::complex<double> ci_tpi = ModuleBase::NEG_IMAG_UNIT * ModuleBase::TWO_PI;
@@ -379,4 +389,4 @@ template <>
 std::complex<double> * Structure_Factor::get_eigts3_data() const
 {
     return this->z_eigts3;
-}
\ No newline at end of file
+}
diff --git a/source/source_pw/module_pwdft/structure_factor.h b/source/source_pw/module_pwdft/structure_factor.h
index 36e9185400..6de0a5b0e3 100644
--- a/source/source_pw/module_pwdft/structure_factor.h
+++ b/source/source_pw/module_pwdft/structure_factor.h
@@ -1,5 +1,5 @@
-#ifndef PLANEWAVE_H
-#define PLANEWAVE_H
+#ifndef STRUCTURE_FACTOR_H 
+#define STRUCTURE_FACTOR_H
 
 #include "source_base/complexmatrix.h"
 #include "source_basis/module_pw/pw_basis_k.h"
@@ -23,9 +23,10 @@ class Structure_Factor
 
 	// structure factor (ntype, ngmc)
     ModuleBase::ComplexMatrix strucFac;
-    void setup_structure_factor(const UnitCell* Ucell,
-                                const Parallel_Grid& pgrid,
-                                const ModulePW::PW_Basis* rho_basis); // Calculate structure factors
+
+	void setup(const UnitCell* Ucell,
+			const Parallel_Grid& pgrid,
+			const ModulePW::PW_Basis* rho_basis); // Calculate structure factors
 
     /// calculate structure factors through Cardinal B-spline interpolation
     void bspline_sf(
@@ -33,6 +34,7 @@ class Structure_Factor
         const UnitCell* Ucell,
         const Parallel_Grid& pgrid,
         const ModulePW::PW_Basis* rho_basis); 
+
     void bsplinecoef(std::complex<double> *b1, std::complex<double> *b2, std::complex<double> *b3, 
                     const int nx, const int ny, const int nz, const int norder);
 
@@ -51,7 +53,9 @@ class Structure_Factor
     // sf with k points
     std::complex<double>* get_sk(const int ik, const int it, const int ia, const ModulePW::PW_Basis_K* wfc_basis) const;
     template <typename FPTYPE, typename Device>
+
     void get_sk(Device* ctx, const int ik, const ModulePW::PW_Basis_K* wfc_basis, std::complex<FPTYPE>* sk) const;
+
     std::complex<double>* get_skq(int ik,
                                   int it,
                                   int ia,
@@ -59,9 +63,16 @@ class Structure_Factor
                                   ModuleBase::Vector3<double> q);
 
   private:
+
     const UnitCell* ucell=nullptr; 
-    std::complex<float> * c_eigts1 = nullptr, * c_eigts2 = nullptr, * c_eigts3 = nullptr;
-    std::complex<double> * z_eigts1 = nullptr, * z_eigts2 = nullptr, * z_eigts3 = nullptr;
+    std::complex<float> * c_eigts1 = nullptr;
+    std::complex<float> * c_eigts2 = nullptr;
+    std::complex<float> * c_eigts3 = nullptr;
+
+    std::complex<double> * z_eigts1 = nullptr;
+    std::complex<double> * z_eigts2 = nullptr;
+    std::complex<double> * z_eigts3 = nullptr;
+
     const ModulePW::PW_Basis* rho_basis = nullptr;
     std::string device = "cpu";
 };
diff --git a/source/source_pw/module_pwdft/test/structure_factor_test.cpp b/source/source_pw/module_pwdft/test/structure_factor_test.cpp
index c0f5078a72..0a1866a3f4 100644
--- a/source/source_pw/module_pwdft/test/structure_factor_test.cpp
+++ b/source/source_pw/module_pwdft/test/structure_factor_test.cpp
@@ -75,7 +75,7 @@ TEST_F(StructureFactorTest, set)
 TEST_F(StructureFactorTest, setup_structure_factor_double)
 {
     rho_basis->npw = 10;
-    SF.setup_structure_factor(ucell,*pgrid,rho_basis);  
+    SF.setup(ucell,*pgrid,rho_basis);  
 
     for (int i=0;i< ucell->nat * (2 * rho_basis->nx + 1);i++) 
     {
@@ -100,7 +100,7 @@ TEST_F(StructureFactorTest, setup_structure_factor_float)
 {
     PARAM.sys.has_float_data = true;
     rho_basis->npw = 10;
-    SF.setup_structure_factor(ucell,*pgrid,rho_basis);  
+    SF.setup(ucell,*pgrid,rho_basis);  
 
     for (int i=0;i< ucell->nat * (2 * rho_basis->nx + 1);i++) 
     {
@@ -125,4 +125,4 @@ int main()
 {
     testing::InitGoogleTest();
     return RUN_ALL_TESTS();
-}
\ No newline at end of file
+}
diff --git a/source/source_relax/test/relax_test.h b/source/source_relax/test/relax_test.h
index d3c6ce5d2f..fa787a909a 100644
--- a/source/source_relax/test/relax_test.h
+++ b/source/source_relax/test/relax_test.h
@@ -26,4 +26,4 @@ void ModuleSymmetry::Symmetry::symmetrize_mat3(ModuleBase::matrix& sigma, const
 void ModuleSymmetry::Symmetry::symmetrize_vec3_nat(double* v)const {};
 Structure_Factor::Structure_Factor() {};
 Structure_Factor::~Structure_Factor(){};
-void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const Parallel_Grid&, const ModulePW::PW_Basis* rho_basis){};
+void Structure_Factor::setup(const UnitCell* Ucell, const Parallel_Grid&, const ModulePW::PW_Basis* rho_basis){};