diff --git a/source/source_estate/module_dm/density_matrix.cpp b/source/source_estate/module_dm/density_matrix.cpp
index 1e59218b0d..e5ad1ba620 100644
--- a/source/source_estate/module_dm/density_matrix.cpp
+++ b/source/source_estate/module_dm/density_matrix.cpp
@@ -5,6 +5,7 @@
 #include "source_base/memory.h"
 #include "source_base/timer.h"
 #include "source_base/tool_title.h"
+#include "source_base/constants.h"
 #include "source_cell/klist.h"
 
 namespace elecstate
@@ -52,22 +53,25 @@ DensityMatrix<TK, TR>::DensityMatrix(const Parallel_Orbitals* paraV_in, const in
 }
 
 
+
 // calculate DMR from DMK using blas for multi-k calculation
-template <>
-void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik_in)
+template <typename TK, typename TR_in, typename TR_out>
+void DensityMatrix_Tools::cal_DMR(
+    const DensityMatrix<TK, TR_in> &dm,
+    std::vector<hamilt::HContainer<TR_out>*> &dmR_out,
+    const int ik_in)
 {
     ModuleBase::TITLE("DensityMatrix", "cal_DMR");
-    using TR = double;
 
     // To check whether DMR has been initialized
-    assert(!this->_DMR.empty() && "DMR has not been initialized!");
+    assert(dmR_out.size()==dm._nspin && "DMR has not been initialized!");
 
     ModuleBase::timer::tick("DensityMatrix", "cal_DMR");
-    const int ld_hk = this->_paraV->nrow;
-    for (int is = 1; is <= this->_nspin; ++is)
+    const int ld_hk = dm._paraV->nrow;
+    for (int is = 1; is <= dm._nspin; ++is)
     {
-        const int ik_begin = this->_nk * (is - 1); // jump this->_nk for spin_down if nspin==2
-        hamilt::HContainer<TR>*const target_DMR = this->_DMR[is - 1];
+        const int ik_begin = dm._nk * (is - 1); // jump dm._nk for spin_down if nspin==2
+        hamilt::HContainer<TR_out>*const target_DMR = dmR_out[is - 1];
         // set zero since this function is called in every scf step
         target_DMR->set_zero();
         #ifdef _OPENMP
@@ -75,25 +79,25 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik_in)
         #endif
         for (int i = 0; i < target_DMR->size_atom_pairs(); ++i)
         {
-            hamilt::AtomPair<TR>& target_ap = target_DMR->get_atom_pair(i);
+            hamilt::AtomPair<TR_out>& target_ap = target_DMR->get_atom_pair(i);
             const int iat1 = target_ap.get_atom_i();
             const int iat2 = target_ap.get_atom_j();
             // get global indexes of whole matrix for each atom in this process
-            const int row_ap = this->_paraV->atom_begin_row[iat1];
-            const int col_ap = this->_paraV->atom_begin_col[iat2];
-            const int row_size = this->_paraV->get_row_size(iat1);
-            const int col_size = this->_paraV->get_col_size(iat2);
+            const int row_ap = dm._paraV->atom_begin_row[iat1];
+            const int col_ap = dm._paraV->atom_begin_col[iat2];
+            const int row_size = dm._paraV->get_row_size(iat1);
+            const int col_size = dm._paraV->get_col_size(iat2);
             const int mat_size = row_size * col_size;
             const int R_size = target_ap.get_R_size();
             assert(row_ap != -1 && col_ap != -1 && "Atom-pair not belong this process");
 
             // calculate kphase and target_mat_ptr
-            std::vector<std::vector<std::complex<double>>> kphase_vec(this->_nk, std::vector<std::complex<double>>(R_size));
-            std::vector<TR*> target_DMR_mat_vec(R_size);
+            std::vector<std::vector<TK>> kphase_vec(dm._nk, std::vector<TK>(R_size));
+            std::vector<TR_out*> target_DMR_mat_vec(R_size);
             for(int iR = 0; iR < R_size; ++iR)
             {
                 const ModuleBase::Vector3<int> R_index = target_ap.get_R_index(iR);
-                hamilt::BaseMatrix<TR>*const target_mat = target_ap.find_matrix(R_index);
+                hamilt::BaseMatrix<TR_out>*const target_mat = target_ap.find_matrix(R_index);
                 #ifdef __DEBUG
                 if (target_mat == nullptr)
                 {
@@ -102,29 +106,28 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik_in)
                 }
                 #endif
                 target_DMR_mat_vec[iR] = target_mat->get_pointer();
-                for(int ik = 0; ik < this->_nk; ++ik)
+                for(int ik = 0; ik < dm._nk; ++ik)
                 {
                     if(ik_in >= 0 && ik_in != ik) { continue; }
                     // cal k_phase
                     // if TK==std::complex<double>, kphase is e^{ikR}
                     const ModuleBase::Vector3<double> dR(R_index[0], R_index[1], R_index[2]);
-                    const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                    const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                     double sinp, cosp;
                     ModuleBase::libm::sincos(arg, &sinp, &cosp);
-                    kphase_vec[ik][iR] = std::complex<double>(cosp, sinp);
+                    kphase_vec[ik][iR] = TK(cosp, sinp);
                 }
             }
 
-            std::vector<std::complex<double>> DMK_mat_trans(mat_size);
-            std::vector<std::complex<double>> tmp_DMR( (PARAM.inp.nspin==4) ? mat_size*R_size : 0);
-            for(int ik = 0; ik < this->_nk; ++ik)
+            std::vector<TK> DMK_mat_trans(mat_size);
+            std::vector<TK> tmp_DMR( (PARAM.inp.nspin==4) ? mat_size*R_size : 0);
+            for(int ik = 0; ik < dm._nk; ++ik)
             {
                 if(ik_in >= 0 && ik_in != ik) { continue; }
-
                 // copy column-major DMK to row-major DMK_mat_trans (for the purpose of computational efficiency)
-                const std::complex<double>*const DMK_mat_ptr
-                    = this->_DMK[ik + ik_begin].data()
-                      + col_ap * this->_paraV->nrow + row_ap;
+                const TK*const DMK_mat_ptr
+                    = dm._DMK[ik + ik_begin].data()
+                      + col_ap * dm._paraV->nrow + row_ap;
                 for(int icol = 0; icol < col_size; ++icol) {
                     for(int irow = 0; irow < row_size; ++irow) {
                         DMK_mat_trans[irow * col_size + icol] = DMK_mat_ptr[icol * ld_hk + irow];
@@ -135,16 +138,10 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik_in)
                 for(int iR = 0; iR < R_size; ++iR)
                 {
                     // (kr+i*ki) * (Dr+i*Di) = (kr*Dr-ki*Di) + i*(kr*Di+ki*Dr)
-                    const std::complex<double> kphase = kphase_vec[ik][iR];
+                    const TK kphase = kphase_vec[ik][iR];
                     if(PARAM.inp.nspin != 4)                // only save real kr*Dr-ki*Di
                     {
-                        TR* target_DMR_mat = target_DMR_mat_vec[iR];
-                        for(int i = 0; i < mat_size; i++)
-                        {
-                            target_DMR_mat[i]
-                                += kphase.real() * DMK_mat_trans[i].real() 
-                                 - kphase.imag() * DMK_mat_trans[i].imag();
-                        }
+                        func_exp_mul_dmk(kphase, DMK_mat_trans, target_DMR_mat_vec[iR]);
                     } else if(PARAM.inp.nspin == 4)
                     {
                         BlasConnector::axpy(mat_size,
@@ -169,11 +166,11 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik_in)
                         step_trace[is * npol + is2] = target_ap.get_col_size() * is + is2;
                 }}
 
-                std::complex<double> tmp[4]{};
+                TK tmp[4]{};
                 for(int iR = 0; iR < R_size; ++iR)
                 {
-                    const std::complex<double>* tmp_DMR_mat = &tmp_DMR[iR * mat_size];
-                    TR* target_DMR_mat = target_DMR_mat_vec[iR];
+                    const TK* tmp_DMR_mat = &tmp_DMR[iR * mat_size];
+                    TR_out* target_DMR_mat = target_DMR_mat_vec[iR];
                     for (int irow = 0; irow < row_size; irow += 2)
                     {
                         for (int icol = 0; icol < col_size; icol += 2)
@@ -183,12 +180,9 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik_in)
                             tmp[1] = tmp_DMR_mat[icol + step_trace[1]];
                             tmp[2] = tmp_DMR_mat[icol + step_trace[2]];
                             tmp[3] = tmp_DMR_mat[icol + step_trace[3]];
-                            // transfer to Pauli matrix and save the real part
-                            // save them back to the target_mat
-                            target_DMR_mat[icol + step_trace[0]] = tmp[0].real() + tmp[3].real();  // (rho_upup + rho_downdown).real()
-                            target_DMR_mat[icol + step_trace[1]] = tmp[1].real() + tmp[2].real();  // (rho_updown + rho_downup).real()
-                            target_DMR_mat[icol + step_trace[2]] = -tmp[1].imag() + tmp[2].imag(); // (i * (rho_updown - rho_downup)).real()
-                            target_DMR_mat[icol + step_trace[3]] = tmp[0].real() - tmp[3].real();  // (rho_upup - rho_downdown).real()
+                            
+                            // transfer to Pauli matrix, save them back to the target_DMR_mat
+                            func_xyz_to_updown(tmp, icol, step_trace, target_DMR_mat);
                         }
                         tmp_DMR_mat += col_size * 2;
                         target_DMR_mat += col_size * 2;
@@ -200,21 +194,39 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik_in)
     ModuleBase::timer::tick("DensityMatrix", "cal_DMR");
 }
 
-// calculate DMR from DMK using blas for multi-k calculation
 template <>
-void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
+void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik_in)
+{
+    DensityMatrix_Tools::cal_DMR(*this, this->_DMR, ik_in);
+}
+
+template <>
+void DensityMatrix<std::complex<double>, std::complex<double>>::cal_DMR(const int ik_in)
+{
+    DensityMatrix_Tools::cal_DMR(*this, this->_DMR, ik_in);
+}
+
+
+
+// calculate DMR from DMK using blas for multi-k calculation
+template <typename TK, typename TR_in, typename TR_out>
+void DensityMatrix_Tools::cal_DMR_td(
+    const DensityMatrix<TK, TR_in> &dm,
+    std::vector<hamilt::HContainer<TR_out>*> &dmR_out,
+    const UnitCell& ucell,
+    const ModuleBase::Vector3<double> At,
+    const int ik_in)
 {
     ModuleBase::TITLE("DensityMatrix", "cal_DMR_td");
-    using TR = double;
     // To check whether DMR has been initialized
-    assert(!this->_DMR.empty() && "DMR has not been initialized!");
+    assert(dmR_out.size()==dm._nspin && "DMR has not been initialized!");
 
     ModuleBase::timer::tick("DensityMatrix", "cal_DMR_td");
-    const int ld_hk = this->_paraV->nrow;
-    for (int is = 1; is <= this->_nspin; ++is)
+    const int ld_hk = dm._paraV->nrow;
+    for (int is = 1; is <= dm._nspin; ++is)
     {
-        const int ik_begin = this->_nk * (is - 1); // jump this->_nk for spin_down if nspin==2
-        hamilt::HContainer<TR>*const target_DMR = this->_DMR[is - 1];
+        const int ik_begin = dm._nk * (is - 1); // jump dm._nk for spin_down if nspin==2
+        hamilt::HContainer<TR_out>*const target_DMR = dmR_out[is - 1];
         // set zero since this function is called in every scf step
         target_DMR->set_zero();
         #ifdef _OPENMP
@@ -222,25 +234,25 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& uce
         #endif
         for (int i = 0; i < target_DMR->size_atom_pairs(); ++i)
         {
-            hamilt::AtomPair<TR>& target_ap = target_DMR->get_atom_pair(i);
+            hamilt::AtomPair<TR_out>& target_ap = target_DMR->get_atom_pair(i);
             const int iat1 = target_ap.get_atom_i();
             const int iat2 = target_ap.get_atom_j();
             // get global indexes of whole matrix for each atom in this process
-            const int row_ap = this->_paraV->atom_begin_row[iat1];
-            const int col_ap = this->_paraV->atom_begin_col[iat2];
-            const int row_size = this->_paraV->get_row_size(iat1);
-            const int col_size = this->_paraV->get_col_size(iat2);
+            const int row_ap = dm._paraV->atom_begin_row[iat1];
+            const int col_ap = dm._paraV->atom_begin_col[iat2];
+            const int row_size = dm._paraV->get_row_size(iat1);
+            const int col_size = dm._paraV->get_col_size(iat2);
             const int mat_size = row_size * col_size;
             const int R_size = target_ap.get_R_size();
             assert(row_ap != -1 && col_ap != -1 && "Atom-pair not belong this process");
 
             // calculate kphase and target_mat_ptr
-            std::vector<std::vector<std::complex<double>>> kphase_vec(this->_nk, std::vector<std::complex<double>>(R_size));
-            std::vector<TR*> target_DMR_mat_vec(R_size);
+            std::vector<std::vector<TK>> kphase_vec(dm._nk, std::vector<TK>(R_size));
+            std::vector<TR_out*> target_DMR_mat_vec(R_size);
             for(int iR = 0; iR < R_size; ++iR)
             {
                 const ModuleBase::Vector3<int> R_index = target_ap.get_R_index(iR);
-                hamilt::BaseMatrix<TR>*const target_mat = target_ap.find_matrix(R_index);
+                hamilt::BaseMatrix<TR_out>*const target_mat = target_ap.find_matrix(R_index);
                 #ifdef __DEBUG
                 if (target_mat == nullptr)
                 {
@@ -249,33 +261,31 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& uce
                 }
                 #endif
                 target_DMR_mat_vec[iR] = target_mat->get_pointer();
-                double arg_td = 0.0;
                 //cal tddft phase for hybrid gauge
-                ModuleBase::Vector3<double> dtau = ucell.cal_dtau(iat1, iat2, R_index);
-                arg_td = At * dtau * ucell.lat0;
-                for(int ik = 0; ik < this->_nk; ++ik)
+                const ModuleBase::Vector3<double> dtau = ucell.cal_dtau(iat1, iat2, R_index);
+                const double arg_td = At * dtau * ucell.lat0;
+                for(int ik = 0; ik < dm._nk; ++ik)
                 {
                     if(ik_in >= 0 && ik_in != ik) { continue; }
                     // cal k_phase
                     // if TK==std::complex<double>, kphase is e^{ikR}
                     const ModuleBase::Vector3<double> dR(R_index[0], R_index[1], R_index[2]);
-                    const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI + arg_td;
+                    const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI + arg_td;
                     double sinp, cosp;
                     ModuleBase::libm::sincos(arg, &sinp, &cosp);
-                    kphase_vec[ik][iR] = std::complex<double>(cosp, sinp);
+                    kphase_vec[ik][iR] = TK(cosp, sinp);
                 }
             }
 
-            std::vector<std::complex<double>> DMK_mat_trans(mat_size);
-            std::vector<std::complex<double>> tmp_DMR( (PARAM.inp.nspin==4) ? mat_size*R_size : 0);
-            for(int ik = 0; ik < this->_nk; ++ik)
+            std::vector<TK> DMK_mat_trans(mat_size);
+            std::vector<TK> tmp_DMR( (PARAM.inp.nspin==4) ? mat_size*R_size : 0);
+            for(int ik = 0; ik < dm._nk; ++ik)
             {
                 if(ik_in >= 0 && ik_in != ik) { continue; }
-
                 // copy column-major DMK to row-major DMK_mat_trans (for the purpose of computational efficiency)
-                const std::complex<double>*const DMK_mat_ptr
-                    = this->_DMK[ik + ik_begin].data()
-                      + col_ap * this->_paraV->nrow + row_ap;
+                const TK*const DMK_mat_ptr
+                    = dm._DMK[ik + ik_begin].data()
+                      + col_ap * dm._paraV->nrow + row_ap;
                 for(int icol = 0; icol < col_size; ++icol) {
                     for(int irow = 0; irow < row_size; ++irow) {
                         DMK_mat_trans[irow * col_size + icol] = DMK_mat_ptr[icol * ld_hk + irow];
@@ -286,16 +296,10 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& uce
                 for(int iR = 0; iR < R_size; ++iR)
                 {
                     // (kr+i*ki) * (Dr+i*Di) = (kr*Dr-ki*Di) + i*(kr*Di+ki*Dr)
-                    const std::complex<double> kphase = kphase_vec[ik][iR];
+                    const TK kphase = kphase_vec[ik][iR];
                     if(PARAM.inp.nspin != 4)                // only save real kr*Dr-ki*Di
                     {
-                        TR* target_DMR_mat = target_DMR_mat_vec[iR];
-                        for(int i = 0; i < mat_size; i++)
-                        {
-                            target_DMR_mat[i]
-                                += kphase.real() * DMK_mat_trans[i].real() 
-                                 - kphase.imag() * DMK_mat_trans[i].imag();
-                        }
+                        func_exp_mul_dmk(kphase, DMK_mat_trans, target_DMR_mat_vec[iR]);
                     } else if(PARAM.inp.nspin == 4)
                     {
                         BlasConnector::axpy(mat_size,
@@ -320,11 +324,11 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& uce
                         step_trace[is * npol + is2] = target_ap.get_col_size() * is + is2;
                 }}
 
-                std::complex<double> tmp[4]{};
+                TK tmp[4]{};
                 for(int iR = 0; iR < R_size; ++iR)
                 {
-                    const std::complex<double>* tmp_DMR_mat = &tmp_DMR[iR * mat_size];
-                    TR* target_DMR_mat = target_DMR_mat_vec[iR];
+                    const TK* tmp_DMR_mat = &tmp_DMR[iR * mat_size];
+                    TR_out* target_DMR_mat = target_DMR_mat_vec[iR];
                     for (int irow = 0; irow < row_size; irow += 2)
                     {
                         for (int icol = 0; icol < col_size; icol += 2)
@@ -334,12 +338,9 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& uce
                             tmp[1] = tmp_DMR_mat[icol + step_trace[1]];
                             tmp[2] = tmp_DMR_mat[icol + step_trace[2]];
                             tmp[3] = tmp_DMR_mat[icol + step_trace[3]];
-                            // transfer to Pauli matrix and save the real part
-                            // save them back to the target_mat
-                            target_DMR_mat[icol + step_trace[0]] = tmp[0].real() + tmp[3].real();  // (rho_upup + rho_downdown).real()
-                            target_DMR_mat[icol + step_trace[1]] = tmp[1].real() + tmp[2].real();  // (rho_updown + rho_downup).real()
-                            target_DMR_mat[icol + step_trace[2]] = -tmp[1].imag() + tmp[2].imag(); // (i * (rho_updown - rho_downup)).real()
-                            target_DMR_mat[icol + step_trace[3]] = tmp[0].real() - tmp[3].real();  // (rho_upup - rho_downdown).real()
+                            
+                            // transfer to Pauli matrix, save them back to the target_DMR_mat
+                            func_xyz_to_updown(tmp, icol, step_trace, target_DMR_mat);
                         }
                         tmp_DMR_mat += col_size * 2;
                         target_DMR_mat += col_size * 2;
@@ -351,17 +352,32 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& uce
     ModuleBase::timer::tick("DensityMatrix", "cal_DMR_td");
 }
 
-// calculate DMR from DMK using blas for multi-k calculation
 template <>
-void DensityMatrix<double, double>::cal_DMR_full(hamilt::HContainer<std::complex<double>>* dmR_out)const{}
+void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
+{
+    DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, ucell, At, ik_in);
+}
+
 template <>
-void DensityMatrix<std::complex<double>, double>::cal_DMR_full(hamilt::HContainer<std::complex<double>>* dmR_out)const
+void DensityMatrix<std::complex<double>, std::complex<double>>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
+{
+    DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, ucell, At, ik_in);
+}
+
+
+
+// calculate DMR from DMK using blas for multi-k calculation
+template <typename TK, typename TR_in, typename TR_out>
+void DensityMatrix_Tools::cal_DMR_full(
+    const DensityMatrix<TK, TR_in> &dm, 
+    hamilt::HContainer<TR_out>* dmR_out,
+    const int ik_in)
 {
     ModuleBase::TITLE("DensityMatrix", "cal_DMR_full");
 
     ModuleBase::timer::tick("DensityMatrix", "cal_DMR_full");
-    const int ld_hk = this->_paraV->nrow;
-    hamilt::HContainer<std::complex<double>>* target_DMR = dmR_out;
+    const int ld_hk = dm._paraV->nrow;
+    hamilt::HContainer<TR_out>* target_DMR = dmR_out;
     // set zero since this function is called in every scf step
     target_DMR->set_zero();
     #ifdef _OPENMP
@@ -369,24 +385,25 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_full(hamilt::HContaine
     #endif
     for (int i = 0; i < target_DMR->size_atom_pairs(); ++i)
     {
-        hamilt::AtomPair<std::complex<double>>& target_ap = target_DMR->get_atom_pair(i);
+        hamilt::AtomPair<TR_out>& target_ap = target_DMR->get_atom_pair(i);
         const int iat1 = target_ap.get_atom_i();
         const int iat2 = target_ap.get_atom_j();
         // get global indexes of whole matrix for each atom in this process
-        const int row_ap = this->_paraV->atom_begin_row[iat1];
-        const int col_ap = this->_paraV->atom_begin_col[iat2];
-        const int row_size = this->_paraV->get_row_size(iat1);
-        const int col_size = this->_paraV->get_col_size(iat2);
+        const int row_ap = dm._paraV->atom_begin_row[iat1];
+        const int col_ap = dm._paraV->atom_begin_col[iat2];
+        const int row_size = dm._paraV->get_row_size(iat1);
+        const int col_size = dm._paraV->get_col_size(iat2);
         const int mat_size = row_size * col_size;
         const int R_size = target_ap.get_R_size();
+        assert(row_ap != -1 && col_ap != -1 && "Atom-pair not belong this process");
 
         // calculate kphase and target_mat_ptr
-        std::vector<std::vector<std::complex<double>>> kphase_vec(this->_nk, std::vector<std::complex<double>>(R_size));
-        std::vector<std::complex<double>*> target_DMR_mat_vec(R_size);
+        std::vector<std::vector<TK>> kphase_vec(dm._nk, std::vector<TK>(R_size));
+        std::vector<TR_out*> target_DMR_mat_vec(R_size);
         for(int iR = 0; iR < R_size; ++iR)
         {
             const ModuleBase::Vector3<int> R_index = target_ap.get_R_index(iR);
-            hamilt::BaseMatrix<std::complex<double>>*const target_mat = target_ap.find_matrix(R_index);
+            hamilt::BaseMatrix<TR_out>*const target_mat = target_ap.find_matrix(R_index);
             #ifdef __DEBUG
             if (target_mat == nullptr)
             {
@@ -395,25 +412,27 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_full(hamilt::HContaine
             }
             #endif
             target_DMR_mat_vec[iR] = target_mat->get_pointer();
-            for(int ik = 0; ik < this->_nk; ++ik)
+            for(int ik = 0; ik < dm._nk; ++ik)
             {
+                if(ik_in >= 0 && ik_in != ik) { continue; }
                 // cal k_phase
                 // if TK==std::complex<double>, kphase is e^{ikR}
                 const ModuleBase::Vector3<double> dR(R_index[0], R_index[1], R_index[2]);
-                const double arg = (this->_kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI;
                 double sinp, cosp;
                 ModuleBase::libm::sincos(arg, &sinp, &cosp);
-                kphase_vec[ik][iR] = std::complex<double>(cosp, sinp);
+                kphase_vec[ik][iR] = TK(cosp, sinp);
             }
         }
 
-        std::vector<std::complex<double>> DMK_mat_trans(mat_size);
-        for(int ik = 0; ik < this->_nk; ++ik)
+        std::vector<TK> DMK_mat_trans(mat_size);
+        for(int ik = 0; ik < dm._nk; ++ik)
         {
+            if(ik_in >= 0 && ik_in != ik) { continue; }
             // copy column-major DMK to row-major DMK_mat_trans (for the purpose of computational efficiency)
-            const std::complex<double>*const DMK_mat_ptr
-                = this->_DMK[ik].data()
-                  + col_ap * this->_paraV->nrow + row_ap;
+            const TK*const DMK_mat_ptr
+                = dm._DMK[ik].data()
+                  + col_ap * dm._paraV->nrow + row_ap;
             for(int icol = 0; icol < col_size; ++icol) {
                 for(int irow = 0; irow < row_size; ++irow) {
                     DMK_mat_trans[irow * col_size + icol] = DMK_mat_ptr[icol * ld_hk + irow];
@@ -421,13 +440,12 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_full(hamilt::HContaine
 
             for(int iR = 0; iR < R_size; ++iR)
             {
-                const std::complex<double> kphase = kphase_vec[ik][iR];
-                std::complex<double>* target_DMR_mat = target_DMR_mat_vec[iR];
+                const TK kphase = kphase_vec[ik][iR];
                 BlasConnector::axpy(mat_size,
                                     kphase,
                                     DMK_mat_trans.data(),
                                     1,
-                                    target_DMR_mat,
+                                    target_DMR_mat_vec[iR],
                                     1);
             }
         }
@@ -435,47 +453,61 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_full(hamilt::HContaine
     ModuleBase::timer::tick("DensityMatrix", "cal_DMR_full");
 }
 
+template <>
+void DensityMatrix<double, double>::cal_DMR_full(
+    hamilt::HContainer<std::complex<double>>* dmR_out,
+    const int ik_in) const{}
+template <>
+void DensityMatrix<std::complex<double>, double>::cal_DMR_full(
+    hamilt::HContainer<std::complex<double>>* dmR_out,
+    const int ik_in) const
+{
+    DensityMatrix_Tools::cal_DMR_full(*this, dmR_out, ik_in);
+}
+
+
+
 // calculate DMR from DMK using blas for gamma-only calculation
 template <>
 void DensityMatrix<double, double>::cal_DMR(const int ik_in)
 {
     ModuleBase::TITLE("DensityMatrix", "cal_DMR");
+    using TK = double;
+    using TR = double;
 
     assert(ik_in == -1 || ik_in == 0);
     assert(this->_nk == 1);
 
     // To check whether DMR has been initialized
-    assert(!this->_DMR.empty() && "DMR has not been initialized!");
+    assert(this->_DMR.size()==this->_nspin && "DMR has not been initialized!");
 
     ModuleBase::timer::tick("DensityMatrix", "cal_DMR");
     const int ld_hk = this->_paraV->nrow;
     for (int is = 1; is <= this->_nspin; ++is)
     {
         const int ik_begin = this->_nk * (is - 1); // jump this->_nk for spin_down if nspin==2
-        hamilt::HContainer<double>*const target_DMR = this->_DMR[is - 1];
+        hamilt::HContainer<TR>*const target_DMR = this->_DMR[is - 1];
         // set zero since this function is called in every scf step
         target_DMR->set_zero();
-
-        // assert(target_DMR->is_gamma_only() == true);
         #ifdef _OPENMP
         #pragma omp parallel for schedule(dynamic)
         #endif
         for (int i = 0; i < target_DMR->size_atom_pairs(); ++i)
         {
-            hamilt::AtomPair<double>& target_ap = target_DMR->get_atom_pair(i);
+            hamilt::AtomPair<TR>& target_ap = target_DMR->get_atom_pair(i);
             const int iat1 = target_ap.get_atom_i();
             const int iat2 = target_ap.get_atom_j();
             // get global indexes of whole matrix for each atom in this process
             const int row_ap = this->_paraV->atom_begin_row[iat1];
             const int col_ap = this->_paraV->atom_begin_col[iat2];
-            assert(row_ap != -1 && col_ap != -1 && "Atom-pair not belong this process");
             const int row_size = this->_paraV->get_row_size(iat1);
             const int col_size = this->_paraV->get_col_size(iat2);
             const int R_size = target_ap.get_R_size();
+            assert(row_ap != -1 && col_ap != -1 && "Atom-pair not belong this process");
             assert(R_size == 1);
             const ModuleBase::Vector3<int> R_index = target_ap.get_R_index(0);
             assert(R_index.x == 0 && R_index.y == 0 && R_index.z == 0);
-            hamilt::BaseMatrix<double>*const target_mat = target_ap.find_matrix(R_index);
+            hamilt::BaseMatrix<TR>*const target_mat = target_ap.find_matrix(R_index);
             #ifdef __DEBUG
             if (target_mat == nullptr)
             {
@@ -484,22 +516,22 @@ void DensityMatrix<double, double>::cal_DMR(const int ik_in)
             }
             #endif
             // k index
-            constexpr double kphase = 1;
+            constexpr TK kphase = 1;
             // transpose DMK col=>row
-            const double* DMK_ptr
+            const TK* DMK_mat_ptr
                 = this->_DMK[0 + ik_begin].data()
                   + col_ap * this->_paraV->nrow + row_ap;
             // set DMR element
-            double* target_DMR_ptr = target_mat->get_pointer();
+            TR* target_DMR_ptr = target_mat->get_pointer();
             for (int mu = 0; mu < row_size; ++mu)
             {
                 BlasConnector::axpy(col_size,
                                     kphase,
-                                    DMK_ptr,
+                                    DMK_mat_ptr,
                                     ld_hk,
                                     target_DMR_ptr,
                                     1);
-                DMK_ptr += 1;
+                DMK_mat_ptr += 1;
                 target_DMR_ptr += col_size;
             }
         }
@@ -507,6 +539,8 @@ void DensityMatrix<double, double>::cal_DMR(const int ik_in)
     ModuleBase::timer::tick("DensityMatrix", "cal_DMR");
 }
 
+
+
 // switch_dmr
 template <typename TK, typename TR>
 void DensityMatrix<TK, TR>::switch_dmr(const int mode)
@@ -584,6 +618,51 @@ void DensityMatrix<TK, TR>::switch_dmr(const int mode)
     }
 }
 
+
+
+template <>
+void DensityMatrix_Tools::func_exp_mul_dmk<double>(const std::complex<double> kphase, const std::vector<std::complex<double>> &DMK_mat_trans, double* target_DMR_mat)
+{
+    const std::size_t mat_size = DMK_mat_trans.size();
+    for(std::size_t i = 0; i < mat_size; i++)
+    {
+        target_DMR_mat[i]
+            += kphase.real() * DMK_mat_trans[i].real() 
+            - kphase.imag() * DMK_mat_trans[i].imag();
+    }
+}
+
+template <>
+void DensityMatrix_Tools::func_exp_mul_dmk<std::complex<double>>(const std::complex<double> kphase, const std::vector<std::complex<double>> &DMK_mat_trans, std::complex<double>* target_DMR_mat)
+{
+    BlasConnector::axpy(DMK_mat_trans.size(),
+                        kphase,
+                        DMK_mat_trans.data(),
+                        1,
+                        target_DMR_mat,
+                        1);
+}
+
+template <>
+void DensityMatrix_Tools::func_xyz_to_updown<double>(const std::complex<double> tmp[4], const int icol, const int step_trace[4], double* target_DMR_mat)
+{
+    target_DMR_mat[icol + step_trace[0]] = tmp[0].real() + tmp[3].real();  // rho_0 = (rho_upup + rho_downdown).real()
+    target_DMR_mat[icol + step_trace[1]] = tmp[1].real() + tmp[2].real();  // rho_x = (rho_updown + rho_downup).real()
+    target_DMR_mat[icol + step_trace[2]] = -tmp[1].imag() + tmp[2].imag(); // rho_y = (i * (rho_updown - rho_downup)).real()
+    target_DMR_mat[icol + step_trace[3]] = tmp[0].real() - tmp[3].real();  // rho_z = (rho_upup - rho_downdown).real()
+}
+
+template <>
+void DensityMatrix_Tools::func_xyz_to_updown<std::complex<double>>(const std::complex<double> tmp[4], const int icol, const int step_trace[4], std::complex<double>* target_DMR_mat)
+{
+    target_DMR_mat[icol + step_trace[0]] = tmp[0] + tmp[3];                                         // rho_0 = (rho_upup + rho_downdown)
+    target_DMR_mat[icol + step_trace[1]] = tmp[1] + tmp[2];                                         // rho_x = (rho_updown + rho_downup)
+    target_DMR_mat[icol + step_trace[2]] = ModuleBase::IMAG_UNIT * (tmp[1].imag() - tmp[2].imag()); // rho_y = (i * (rho_updown - rho_downup))
+    target_DMR_mat[icol + step_trace[3]] = tmp[0] - tmp[3];                                         // rho_z = (rho_upup - rho_downdown)
+}
+
+
+
 // T of HContainer can be double or complex<double>
 template class DensityMatrix<double, double>;               // Gamma-Only case
 template class DensityMatrix<std::complex<double>, double>; // Multi-k case
diff --git a/source/source_estate/module_dm/density_matrix.h b/source/source_estate/module_dm/density_matrix.h
index 605e7b6acd..a8b0e1c4ec 100644
--- a/source/source_estate/module_dm/density_matrix.h
+++ b/source/source_estate/module_dm/density_matrix.h
@@ -31,6 +31,40 @@ struct ShiftRealComplex<std::complex<double>>
 	using type = double;
 };
 
+
+    template <typename TK, typename TR> class DensityMatrix;
+
+// DensityMatrix<complex<double>,TR>::cal_DMR() is illegal in C++, so DensityMatrix_Tools is used instead.
+namespace DensityMatrix_Tools
+{
+    template <typename TK, typename TR_in, typename TR_out>
+    extern void cal_DMR(
+        const DensityMatrix<TK, TR_in> &dm,
+        std::vector<hamilt::HContainer<TR_out>*> &dmR_out,
+        const int ik_in);
+
+    template <typename TK, typename TR_in, typename TR_out>
+    extern void cal_DMR_td(
+        const DensityMatrix<TK, TR_in> &dm,
+        std::vector<hamilt::HContainer<TR_out>*> &dmR_out,
+        const UnitCell& ucell,
+        const ModuleBase::Vector3<double> At,
+        const int ik_in);
+
+    template <typename TK, typename TR_in, typename TR_out>
+    extern void cal_DMR_full(
+        const DensityMatrix<TK, TR_in> &dm, 
+        hamilt::HContainer<TR_out>* dmR_out,
+        const int ik_in);
+
+    template <typename TR>
+    extern void func_exp_mul_dmk(const std::complex<double> kphase, const std::vector<std::complex<double>> &DMK_mat_trans, TR* target_DMR_mat);
+
+    template <typename TR>
+    extern void func_xyz_to_updown(const std::complex<double> tmp[4], const int icol, const int step_trace[4], TR* target_DMR_mat);
+}
+
+
 template <typename TK, typename TR>
 class DensityMatrix
 {
@@ -179,6 +213,7 @@ class DensityMatrix
 
     /**
      * @brief calculate density matrix DMR from dm(k) using blas::axpy
+     * @param ik_in
      * if ik_in < 0, calculate all k-points
      * if ik_in >= 0, calculate only one k-point without summing over k-points
      */
@@ -186,6 +221,7 @@ class DensityMatrix
 
     /**
      * @brief calculate density matrix DMR with additional vector potential phase, used for hybrid gauge tddft
+     * @param ik_in
      * if ik_in < 0, calculate all k-points
      * if ik_in >= 0, calculate only one k-point without summing over k-points
      */
@@ -195,8 +231,11 @@ class DensityMatrix
      * @brief calculate complex density matrix DMR with both real and imaginary part for noncollinear-spin calculation
      * the stored dm(k) has been used to calculate the passin DMR
      * @param dmR_out pointer of HContainer object to store the calculated complex DMR
+     * @param ik_in
+     * if ik_in < 0, calculate all k-points
+     * if ik_in >= 0, calculate only one k-point without summing over k-points
      */
-    void cal_DMR_full(hamilt::HContainer<std::complex<double>>* dmR_out) const;
+    void cal_DMR_full(hamilt::HContainer<std::complex<double>>* dmR_out, const int ik_in = -1) const;
 
     /**
      * @brief (Only nspin=2) switch DMR to total density matrix or magnetization density matrix
@@ -287,6 +326,9 @@ class DensityMatrix
     std::vector<TR> dmr_origin_;
     TR* dmr_tmp_ = nullptr;
 
+    friend void DensityMatrix_Tools::cal_DMR<TK,TR>(const DensityMatrix<TK, TR> &dm, std::vector<hamilt::HContainer<TR>*> &dmR_out, const int ik_in);
+    friend void DensityMatrix_Tools::cal_DMR_td<TK,TR>(const DensityMatrix<TK, TR> &dm, std::vector<hamilt::HContainer<TR>*> &dmR_out, const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in);
+    friend void DensityMatrix_Tools::cal_DMR_full<TK,TR>(const DensityMatrix<TK, TR> &dm, hamilt::HContainer<std::complex<double>>* dmR_out, const int ik_in);
 };
 
 } // namespace elecstate
diff --git a/source/source_lcao/module_lr/CMakeLists.txt b/source/source_lcao/module_lr/CMakeLists.txt
index 3278252c84..331e7dc22e 100644
--- a/source/source_lcao/module_lr/CMakeLists.txt
+++ b/source/source_lcao/module_lr/CMakeLists.txt
@@ -11,7 +11,6 @@ if(ENABLE_LCAO)
     ao_to_mo_transformer/ao_to_mo_serial.cpp
     dm_trans/dm_trans_parallel.cpp
     dm_trans/dm_trans_serial.cpp
-    dm_trans/dmr_complex.cpp
     operator_casida/operator_lr_hxc.cpp
     operator_casida/operator_lr_exx.cpp
     potentials/pot_hxc_lrtd.cpp