diff --git a/source/source_estate/elecstate_lcao.cpp b/source/source_estate/elecstate_lcao.cpp
index 14562fca87..7bd32542ae 100644
--- a/source/source_estate/elecstate_lcao.cpp
+++ b/source/source_estate/elecstate_lcao.cpp
@@ -1,6 +1,5 @@
-#include "elecstate_lcao.h"
-
-#include "cal_dm.h"
+#include "source_estate/elecstate_lcao.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"
@@ -31,11 +30,12 @@ double ElecStateLCAO<std::complex<double>>::get_spin_constrain_energy()
     return sc.cal_escon();
 }
 
-#ifdef __PEXSI
 template <>
-void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM, std::vector<double*> pexsi_EDM)
+void ElecStateLCAO<double>::dm2rho(std::vector<double*> pexsi_DM, 
+		std::vector<double*> pexsi_EDM,
+		DensityMatrix<double, double>* dm)
 {
-    ModuleBase::timer::tick("ElecStateLCAO", "dm2Rho");
+    ModuleBase::timer::tick("ElecStateLCAO", "dm2rho");
 
     int nspin = PARAM.inp.nspin;
     if (PARAM.inp.nspin == 4)
@@ -43,13 +43,15 @@ void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM, std::vector<do
         nspin = 1;
     }
 
-    this->get_DM()->pexsi_EDM = pexsi_EDM;
+#ifdef __PEXSI
+    dm->pexsi_EDM = pexsi_EDM;
+#endif
 
     for (int is = 0; is < nspin; is++)
     {
-        this->DM->set_DMK_pointer(is, pexsi_DM[is]);
+        dm->set_DMK_pointer(is, pexsi_DM[is]);
     }
-    DM->cal_DMR();
+    dm->cal_DMR();
 
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
@@ -58,30 +60,30 @@ void ElecStateLCAO<double>::dm2Rho(std::vector<double*> pexsi_DM, std::vector<do
     }
 
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
-    ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho);
+    ModuleGint::cal_gint_rho(dm->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho);
     if (XC_Functional::get_ked_flag())
     {
         for (int is = 0; is < PARAM.inp.nspin; is++)
         {
             ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0], this->charge->nrxx);
         }
-        ModuleGint::cal_gint_tau(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r);
+        ModuleGint::cal_gint_tau(dm->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r);
     }
 
     this->charge->renormalize_rho();
 
-    ModuleBase::timer::tick("ElecStateLCAO", "dm2Rho");
+    ModuleBase::timer::tick("ElecStateLCAO", "dm2rho");
     return;
 }
 
 template <>
 void ElecStateLCAO<std::complex<double>>::dm2rho(std::vector<std::complex<double>*> pexsi_DM,
-                                                  std::vector<std::complex<double>*> pexsi_EDM)
+		std::vector<std::complex<double>*> pexsi_EDM,
+		DensityMatrix<std::complex<double>, double>* dm)
 {
     ModuleBase::WARNING_QUIT("ElecStateLCAO", "pexsi is not completed for multi-k case");
 }
 
-#endif
 
 template class ElecStateLCAO<double>;               // Gamma_only case
 template class ElecStateLCAO<std::complex<double>>; // multi-k case
diff --git a/source/source_estate/elecstate_lcao.h b/source/source_estate/elecstate_lcao.h
index 5d2d9ac7f6..3cf06641e3 100644
--- a/source/source_estate/elecstate_lcao.h
+++ b/source/source_estate/elecstate_lcao.h
@@ -37,7 +37,6 @@ class ElecStateLCAO : public ElecState
 
     double get_spin_constrain_energy() override;
 
-#ifdef __PEXSI
     // use for pexsi
 
     /**
@@ -46,8 +45,9 @@ class ElecStateLCAO : public ElecState
      * @param pexsi_EDM: pointers of energy-weighed density matrix (EDMK) calculated by pexsi, needed by MD, will be
      * stored in DensityMatrix::pexsi_EDM
      */
-    void dm2rho(std::vector<TK*> pexsi_DM, std::vector<TK*> pexsi_EDM);
-#endif
+	void dm2rho(std::vector<TK*> pexsi_DM, 
+			std::vector<TK*> pexsi_EDM, 
+			DensityMatrix<TK, double>* dm);
 
 };
 
diff --git a/source/source_hsolver/hsolver_lcao.cpp b/source/source_hsolver/hsolver_lcao.cpp
index 87eec7475e..1c4c5a5d61 100644
--- a/source/source_hsolver/hsolver_lcao.cpp
+++ b/source/source_hsolver/hsolver_lcao.cpp
@@ -112,7 +112,7 @@ void HSolverLCAO<TK, Device>::solve(hamilt::Hamilt<TK>* pHamilt,
     else if (this->method == "pexsi")
     {
 #ifdef __PEXSI // other purification methods should follow this routine
-        DiagoPexsi<T> pe(ParaV);
+        DiagoPexsi<TK> pe(ParaV);
         for (int ik = 0; ik < psi.get_nk(); ++ik)
         {
             /// update H(k) for each k point
@@ -121,10 +121,10 @@ void HSolverLCAO<TK, Device>::solve(hamilt::Hamilt<TK>* pHamilt,
             // solve eigenvector and eigenvalue for H(k)
             pe.diag(pHamilt, psi, nullptr);
         }
-        auto _pes = dynamic_cast<elecstate::ElecStateLCAO<T>*>(pes);
+        auto _pes = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(pes);
         pes->f_en.eband = pe.totalFreeEnergy;
         // maybe eferm could be dealt with in the future
-        _pes->dm2rho(dm, pe.EDM);
+        _pes->dm2rho(pe.DM, pe.EDM, &dm);
 #endif
     }
 
diff --git a/source/source_lcao/edm.cpp b/source/source_lcao/edm.cpp
index 2a02afcaaf..a6ea90b6c7 100644
--- a/source/source_lcao/edm.cpp
+++ b/source/source_lcao/edm.cpp
@@ -31,10 +31,10 @@ elecstate::DensityMatrix<double, double> Force_LCAO<double>::cal_edm(const elecs
 #ifdef __PEXSI
     if (PARAM.inp.ks_solver == "pexsi")
     {
-        auto pes = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec);
+        // auto pes = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec);
         for (int ik = 0; ik < nspin; ik++)
         {
-            edm.set_DMK_pointer(ik, pes->get_DM()->pexsi_EDM[ik]);
+            edm.set_DMK_pointer(ik, dm.pexsi_EDM[ik]);
         }
         
     }
@@ -47,7 +47,8 @@ elecstate::DensityMatrix<double, double> Force_LCAO<double>::cal_edm(const elecs
 }
 
 template<>
-elecstate::DensityMatrix<std::complex<double>, double> Force_LCAO<std::complex<double>>::cal_edm(const elecstate::ElecState* pelec,
+elecstate::DensityMatrix<std::complex<double>, double> Force_LCAO<std::complex<double>>::cal_edm(
+    const elecstate::ElecState* pelec,
     const psi::Psi<std::complex<double>>& psi,
     const elecstate::DensityMatrix<std::complex<double>, double>& dm,
     const K_Vectors& kv,