diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
index c00c27837f..01b98ce492 100644
--- a/source/module_esolver/esolver_ks_lcao.cpp
+++ b/source/module_esolver/esolver_ks_lcao.cpp
@@ -1039,14 +1039,14 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
 #ifdef __EXX
     if (GlobalC::exx_info.info_global.cal_exx) // Peize Lin add if 2022.11.14
     {
-        const std::string file_name_exx = GlobalV::global_out_dir + "HexxR_" + std::to_string(GlobalV::MY_RANK);
+        const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK);
         if (GlobalC::exx_info.info_ri.real_number)
         {
-            this->exd->write_Hexxs(file_name_exx);
+            this->exd->write_Hexxs_csr(file_name_exx, GlobalC::ucell);
         }
         else
         {
-            this->exc->write_Hexxs(file_name_exx);
+            this->exc->write_Hexxs_csr(file_name_exx, GlobalC::ucell);
         }
     }
 #endif
diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp
index 290feed8b3..09b6183ffb 100644
--- a/source/module_esolver/esolver_ks_lcao_elec.cpp
+++ b/source/module_esolver/esolver_ks_lcao_elec.cpp
@@ -529,11 +529,11 @@ void ESolver_KS_LCAO<TK, TR>::nscf()
     if (GlobalC::exx_info.info_global.cal_exx)
     {
         // GlobalC::exx_lcao.cal_exx_elec_nscf(this->LOWF.ParaV[0]);
-        const std::string file_name_exx = GlobalV::global_out_dir + "HexxR_" + std::to_string(GlobalV::MY_RANK);
+        const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK);
         if (GlobalC::exx_info.info_ri.real_number)
-            this->exd->read_Hexxs(file_name_exx);
+            this->exd->read_Hexxs_csr(file_name_exx, GlobalC::ucell);
         else
-            this->exc->read_Hexxs(file_name_exx);
+            this->exc->read_Hexxs_csr(file_name_exx, GlobalC::ucell);
 
         hamilt::HamiltLCAO<TK, TR>* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
         auto exx = new hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>>(&this->LM,
diff --git a/source/module_io/csr_reader.cpp b/source/module_io/csr_reader.cpp
index 5d6ee372b5..9c5ce359c1 100644
--- a/source/module_io/csr_reader.cpp
+++ b/source/module_io/csr_reader.cpp
@@ -135,6 +135,6 @@ int csrFileReader<T>::getStep() const
 // T of AtomPair can be double
 template class csrFileReader<double>;
 // ToDo: T of AtomPair can be std::complex<double>
-// template class csrFileReader<std::complex<double>>;
+template class csrFileReader<std::complex<double>>;
 
 } // namespace ModuleIO
\ No newline at end of file
diff --git a/source/module_io/single_R_io.cpp b/source/module_io/single_R_io.cpp
index c278d17782..d802cd8371 100644
--- a/source/module_io/single_R_io.cpp
+++ b/source/module_io/single_R_io.cpp
@@ -3,132 +3,35 @@
 #include "module_base/global_function.h"
 #include "module_base/global_variable.h"
 
-void ModuleIO::output_single_R(std::ofstream &ofs, const std::map<size_t, std::map<size_t, double>> &XR, const double &sparse_threshold, const bool &binary, const Parallel_Orbitals &pv)
+inline void write_data(std::ofstream& ofs, const double& data)
 {
-    double *line = nullptr;
-    std::vector<int> indptr;
-    indptr.reserve(GlobalV::NLOCAL + 1);
-    indptr.push_back(0);
-
-    std::stringstream tem1;
-    tem1 << GlobalV::global_out_dir << "temp_sparse_indices.dat";
-    std::ofstream ofs_tem1;
-    std::ifstream ifs_tem1;
-
-    if (GlobalV::DRANK == 0)
-    {
-        if (binary)
-        {
-            ofs_tem1.open(tem1.str().c_str(), std::ios::binary);
-        }
-        else
-        {
-            ofs_tem1.open(tem1.str().c_str());
-        }
-    }
-
-    line = new double[GlobalV::NLOCAL];
-    for(int row = 0; row < GlobalV::NLOCAL; ++row)
-    {
-        // line = new double[GlobalV::NLOCAL];
-        ModuleBase::GlobalFunc::ZEROS(line, GlobalV::NLOCAL);
-
-        if(pv.global2local_row(row) >= 0)
-        {
-            auto iter = XR.find(row);
-            if (iter != XR.end())
-            {
-                for (auto &value : iter->second)
-                {
-                    line[value.first] = value.second;
-                }
-            }
-        }
-
-        Parallel_Reduce::reduce_all(line, GlobalV::NLOCAL);
-
-        if(GlobalV::DRANK == 0)
-        {
-            int nonzeros_count = 0;
-            for (int col = 0; col < GlobalV::NLOCAL; ++col)
-            {
-                if (std::abs(line[col]) > sparse_threshold)
-                {
-                    if (binary)
-                    {
-                        ofs.write(reinterpret_cast<char *>(&line[col]), sizeof(double));
-                        ofs_tem1.write(reinterpret_cast<char *>(&col), sizeof(int));
-                    }
-                    else
-                    {
-                        ofs << " " << std::fixed << std::scientific << std::setprecision(8) << line[col];
-                        ofs_tem1 << " " << col;
-                    }
-
-                    nonzeros_count++;
-
-                }
-
-            }
-            nonzeros_count += indptr.back();
-            indptr.push_back(nonzeros_count);
-        }
-
-        // delete[] line;
-        // line = nullptr;
-
-    }
-
-    delete[] line;
-    line = nullptr;
-
-    if (GlobalV::DRANK == 0)
-    {
-        if (binary)
-        {
-            ofs_tem1.close();
-            ifs_tem1.open(tem1.str().c_str(), std::ios::binary);
-            ofs << ifs_tem1.rdbuf();
-            ifs_tem1.close();
-            for (auto &i : indptr)
-            {
-                ofs.write(reinterpret_cast<char *>(&i), sizeof(int));
-            }
-        }
-        else
-        {
-            ofs << std::endl;
-            ofs_tem1 << std::endl;
-            ofs_tem1.close();
-            ifs_tem1.open(tem1.str().c_str());
-            ofs << ifs_tem1.rdbuf();
-            ifs_tem1.close();
-            for (auto &i : indptr)
-            {
-                ofs << " " << i;
-            }
-            ofs << std::endl;
-        }
-
-        std::remove(tem1.str().c_str());
-
-    }
-
+    ofs << " " << std::fixed << std::scientific << std::setprecision(8) << data;
+}
+inline void write_data(std::ofstream& ofs, const std::complex<double>& data)
+{
+    ofs << " (" << std::fixed << std::scientific << std::setprecision(8) << data.real() << ","
+        << std::fixed << std::scientific << std::setprecision(8) << data.imag() << ")";
 }
 
-void ModuleIO::output_soc_single_R(std::ofstream &ofs, const std::map<size_t, std::map<size_t, std::complex<double>>> &XR, const double &sparse_threshold, const bool &binary, const Parallel_Orbitals &pv)
+template<typename T>
+void ModuleIO::output_single_R(std::ofstream& ofs,
+    const std::map<size_t, std::map<size_t, T>>& XR,
+    const double& sparse_threshold,
+    const bool& binary,
+    const Parallel_Orbitals& pv,
+    const bool& reduce)
 {
-    std::complex<double> *line = nullptr;
+    T* line = nullptr;
     std::vector<int> indptr;
     indptr.reserve(GlobalV::NLOCAL + 1);
     indptr.push_back(0);
 
     std::stringstream tem1;
-    tem1 << GlobalV::global_out_dir << "temp_sparse_indices.dat";
+    tem1 << GlobalV::global_out_dir << std::to_string(GlobalV::DRANK) + "temp_sparse_indices.dat";
     std::ofstream ofs_tem1;
     std::ifstream ifs_tem1;
 
-    if (GlobalV::DRANK == 0)
+    if (!reduce || GlobalV::DRANK == 0)
     {
         if (binary)
         {
@@ -140,13 +43,12 @@ void ModuleIO::output_soc_single_R(std::ofstream &ofs, const std::map<size_t, st
         }
     }
 
-    line = new std::complex<double>[GlobalV::NLOCAL];
+    line = new T[GlobalV::NLOCAL];
     for(int row = 0; row < GlobalV::NLOCAL; ++row)
     {
-        // line = new std::complex<double>[GlobalV::NLOCAL];
         ModuleBase::GlobalFunc::ZEROS(line, GlobalV::NLOCAL);
 
-        if(pv.global2local_row(row) >= 0)
+        if (!reduce || pv.global2local_row(row) >= 0)
         {
             auto iter = XR.find(row);
             if (iter != XR.end())
@@ -158,9 +60,9 @@ void ModuleIO::output_soc_single_R(std::ofstream &ofs, const std::map<size_t, st
             }
         }
 
-        Parallel_Reduce::reduce_all(line, GlobalV::NLOCAL);
+        if (reduce)Parallel_Reduce::reduce_all(line, GlobalV::NLOCAL);
 
-        if (GlobalV::DRANK == 0)
+        if (!reduce || GlobalV::DRANK == 0)
         {
             int nonzeros_count = 0;
             for (int col = 0; col < GlobalV::NLOCAL; ++col)
@@ -169,13 +71,12 @@ void ModuleIO::output_soc_single_R(std::ofstream &ofs, const std::map<size_t, st
                 {
                     if (binary)
                     {
-                        ofs.write(reinterpret_cast<char *>(&line[col]), sizeof(std::complex<double>));
+                        ofs.write(reinterpret_cast<char*>(&line[col]), sizeof(T));
                         ofs_tem1.write(reinterpret_cast<char *>(&col), sizeof(int));
                     }
                     else
                     {
-                        ofs << " (" << std::fixed << std::scientific << std::setprecision(8) << line[col].real() << "," 
-                                    << std::fixed << std::scientific << std::setprecision(8) << line[col].imag() << ")";
+                        write_data(ofs, line[col]);
                         ofs_tem1 << " " << col;
                     }
 
@@ -196,7 +97,7 @@ void ModuleIO::output_soc_single_R(std::ofstream &ofs, const std::map<size_t, st
     delete[] line;
     line = nullptr;
 
-    if (GlobalV::DRANK == 0)
+    if (!reduce || GlobalV::DRANK == 0)
     {
         if (binary)
         {
@@ -226,5 +127,17 @@ void ModuleIO::output_soc_single_R(std::ofstream &ofs, const std::map<size_t, st
 
         std::remove(tem1.str().c_str());
     }
-
 }
+
+template void ModuleIO::output_single_R<double>(std::ofstream& ofs,
+    const std::map<size_t, std::map<size_t, double>>& XR,
+    const double& sparse_threshold,
+    const bool& binary,
+    const Parallel_Orbitals& pv,
+    const bool& reduce);
+template void ModuleIO::output_single_R<std::complex<double>>(std::ofstream& ofs,
+    const std::map<size_t, std::map<size_t, std::complex<double>>>& XR,
+    const double& sparse_threshold,
+    const bool& binary,
+    const Parallel_Orbitals& pv,
+    const bool& reduce);
\ No newline at end of file
diff --git a/source/module_io/single_R_io.h b/source/module_io/single_R_io.h
index d45539714a..76abc590ad 100644
--- a/source/module_io/single_R_io.h
+++ b/source/module_io/single_R_io.h
@@ -5,8 +5,13 @@
 
 namespace ModuleIO
 {
-	void output_single_R(std::ofstream &ofs, const std::map<size_t, std::map<size_t, double>> &XR, const double &sparse_threshold, const bool &binary, const Parallel_Orbitals &pv);
-    	void output_soc_single_R(std::ofstream &ofs, const std::map<size_t, std::map<size_t, std::complex<double>>> &XR, const double &sparse_threshold, const bool &binary, const Parallel_Orbitals &pv);
+    template <typename T>
+    void output_single_R(std::ofstream& ofs,
+        const std::map<size_t, std::map<size_t, T>>& XR,
+        const double& sparse_threshold,
+        const bool& binary,
+        const Parallel_Orbitals& pv,
+        const bool& reduce = true);
 }
 
 #endif
diff --git a/source/module_io/sparse_matrix.cpp b/source/module_io/sparse_matrix.cpp
index b7585bd478..54a564ef98 100644
--- a/source/module_io/sparse_matrix.cpp
+++ b/source/module_io/sparse_matrix.cpp
@@ -92,7 +92,7 @@ void SparseMatrix<T>::readCSR(const std::vector<T>& values,
 
 // define the operator to index a matrix element
 template <typename T>
-T SparseMatrix<T>::operator()(int row, int col)
+T SparseMatrix<T>::operator()(int row, int col) const
 {
     if (row < 0 || row >= _rows || col < 0 || col >= _cols)
     {
diff --git a/source/module_io/sparse_matrix.h b/source/module_io/sparse_matrix.h
index 63b7f6343a..79e427a394 100644
--- a/source/module_io/sparse_matrix.h
+++ b/source/module_io/sparse_matrix.h
@@ -66,7 +66,7 @@ class SparseMatrix
     }
 
     // define the operator to index a matrix element
-    T operator()(int row, int col);
+    T operator()(int row, int col)const;
 
     // set the threshold
     void setSparseThreshold(double sparse_threshold)
diff --git a/source/module_io/write_HS_R.cpp b/source/module_io/write_HS_R.cpp
index 76bf11c7ca..5869b6ee5f 100644
--- a/source/module_io/write_HS_R.cpp
+++ b/source/module_io/write_HS_R.cpp
@@ -125,7 +125,7 @@ void ModuleIO::output_S_R(
     ModuleBase::timer::tick("ModuleIO","output_S_R"); 
 
     UHM.cal_SR_sparse(sparse_threshold, p_ham);
-    ModuleIO::save_SR_sparse(*UHM.LM, sparse_threshold, binary, SR_filename);
+    ModuleIO::save_sparse(UHM.LM->SR_sparse, UHM.LM->all_R_coor, sparse_threshold, binary, SR_filename, *UHM.LM->ParaV, "S", 0);
     UHM.destroy_all_HSR_sparse();
 
     ModuleBase::timer::tick("ModuleIO","output_S_R");
@@ -154,7 +154,7 @@ void ModuleIO::output_T_R(
     }
 
     UHM.cal_TR_sparse(sparse_threshold);
-    ModuleIO::save_TR_sparse(istep, *UHM.LM, sparse_threshold, binary, sst.str().c_str());
+    ModuleIO::save_sparse(UHM.LM->TR_sparse, UHM.LM->all_R_coor, sparse_threshold, binary, sst.str().c_str(), *UHM.LM->ParaV, "T", istep);
     UHM.destroy_TR_sparse();
 
     ModuleBase::timer::tick("ModuleIO","output_T_R");
diff --git a/source/module_io/write_HS_sparse.cpp b/source/module_io/write_HS_sparse.cpp
index 89a10e1dfb..4abd1b6acf 100644
--- a/source/module_io/write_HS_sparse.cpp
+++ b/source/module_io/write_HS_sparse.cpp
@@ -283,7 +283,7 @@ void ModuleIO::save_HSR_sparse(
                 }
                 else
                 {
-                    output_soc_single_R(g1[ispin], HR_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
+                    output_single_R(g1[ispin], HR_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
                 }
             }
         }
@@ -312,7 +312,7 @@ void ModuleIO::save_HSR_sparse(
             }
             else
             {
-                output_soc_single_R(g2, SR_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
+                output_single_R(g2, SR_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
             }
         }
 
@@ -338,292 +338,6 @@ void ModuleIO::save_HSR_sparse(
     return;
 }
 
-void ModuleIO::save_SR_sparse(
-    LCAO_Matrix &lm,
-    const double& sparse_threshold,
-    const bool &binary,  
-    const std::string &SR_filename
-)
-{
-    ModuleBase::TITLE("ModuleIO","save_SR_sparse");
-    ModuleBase::timer::tick("ModuleIO","save_SR_sparse");
-
-    auto &all_R_coor_ptr = lm.all_R_coor;
-    auto &SR_sparse_ptr = lm.SR_sparse;
-    auto &SR_soc_sparse_ptr = lm.SR_soc_sparse;
-
-    int total_R_num = all_R_coor_ptr.size();
-    int output_R_number = 0;
-    int *S_nonzero_num = nullptr;
-
-    S_nonzero_num = new int[total_R_num];
-    ModuleBase::GlobalFunc::ZEROS(S_nonzero_num, total_R_num);
-
-    int count = 0;
-    for (auto &R_coor : all_R_coor_ptr)
-    {
-        if (GlobalV::NSPIN != 4)
-        {
-            auto iter = SR_sparse_ptr.find(R_coor);
-            if (iter != SR_sparse_ptr.end())
-            {
-                for (auto &row_loop : iter->second)
-                {
-                    S_nonzero_num[count] += row_loop.second.size();
-                }
-            }
-        }
-        else
-        {
-            auto iter = SR_soc_sparse_ptr.find(R_coor);
-            if (iter != SR_soc_sparse_ptr.end())
-            {
-                for (auto &row_loop : iter->second)
-                {
-                    S_nonzero_num[count] += row_loop.second.size();
-                }
-            }
-        }
-
-        count++;
-    }
-
-    Parallel_Reduce::reduce_all(S_nonzero_num, total_R_num);
-
-    for (int index = 0; index < total_R_num; ++index)
-    {
-        if (S_nonzero_num[index] != 0)
-        {
-            output_R_number++;
-        }
-    }
-
-    std::stringstream sss;
-    sss << SR_filename;
-    std::ofstream g2;
-
-    if(GlobalV::DRANK==0)
-    {
-        if (binary)
-        {
-            g2.open(sss.str().c_str(), std::ios::binary);
-            g2.write(reinterpret_cast<char *>(0), sizeof(int));
-            g2.write(reinterpret_cast<char *>(&GlobalV::NLOCAL), sizeof(int));
-            g2.write(reinterpret_cast<char *>(&output_R_number), sizeof(int));
-        }
-        else
-        {
-            g2.open(sss.str().c_str());
-            g2 << "STEP: " << 0 << std::endl;
-            g2 << "Matrix Dimension of S(R): " << GlobalV::NLOCAL <<std::endl;
-            g2 << "Matrix number of S(R): " << output_R_number << std::endl;
-        }
-    }
-
-    count = 0;
-    for (auto &R_coor : all_R_coor_ptr)
-    {
-        int dRx = R_coor.x;
-        int dRy = R_coor.y;
-        int dRz = R_coor.z;
-
-        if (S_nonzero_num[count] == 0)
-        {
-            count++;
-            continue;
-        }
-
-        if (GlobalV::DRANK == 0)
-        {
-            if (binary)
-            {
-                g2.write(reinterpret_cast<char *>(&dRx), sizeof(int));
-                g2.write(reinterpret_cast<char *>(&dRy), sizeof(int));
-                g2.write(reinterpret_cast<char *>(&dRz), sizeof(int));
-                g2.write(reinterpret_cast<char *>(&S_nonzero_num[count]), sizeof(int));
-            }
-            else
-            {
-                g2 << dRx << " " << dRy << " " << dRz << " " << S_nonzero_num[count] << std::endl;
-            }
-        }
-
-        if (GlobalV::NSPIN != 4)
-        {
-            output_single_R(g2, SR_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
-        }
-        else
-        {
-            output_soc_single_R(g2, SR_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
-        }
-
-        count++;
-
-    }
-
-    if(GlobalV::DRANK==0) 
-    {
-        g2.close();
-    }
-
-    delete[] S_nonzero_num;
-    S_nonzero_num = nullptr;
-
-    ModuleBase::timer::tick("ModuleIO","save_SR_sparse");
-    return;
-}
-
-void ModuleIO::save_TR_sparse(
-    const int &istep,
-    LCAO_Matrix &lm,
-    const double& sparse_threshold,
-    const bool &binary,  
-    const std::string &TR_filename
-)
-{
-    ModuleBase::TITLE("ModuleIO","save_TR_sparse");
-    ModuleBase::timer::tick("ModuleIO","save_TR_sparse");
-
-    auto &all_R_coor_ptr = lm.all_R_coor;
-    auto &TR_sparse_ptr = lm.TR_sparse;
-    auto &TR_soc_sparse_ptr = lm.TR_soc_sparse;
-
-    int total_R_num = all_R_coor_ptr.size();
-    int output_R_number = 0;
-    int *T_nonzero_num = nullptr;
-    int step = istep;
-
-    T_nonzero_num = new int[total_R_num];
-    ModuleBase::GlobalFunc::ZEROS(T_nonzero_num, total_R_num);
-
-    int count = 0;
-    for (auto &R_coor : all_R_coor_ptr)
-    {
-        if (GlobalV::NSPIN != 4)
-        {
-            auto iter = TR_sparse_ptr.find(R_coor);
-            if (iter != TR_sparse_ptr.end())
-            {
-                for (auto &row_loop : iter->second)
-                {
-                    T_nonzero_num[count] += row_loop.second.size();
-                }
-            }
-        }
-        else
-        {
-            auto iter = TR_soc_sparse_ptr.find(R_coor);
-            if (iter != TR_soc_sparse_ptr.end())
-            {
-                for (auto &row_loop : iter->second)
-                {
-                    T_nonzero_num[count] += row_loop.second.size();
-                }
-            }
-        }
-
-        count++;
-    }
-
-    Parallel_Reduce::reduce_all(T_nonzero_num, total_R_num);
-
-    for (int index = 0; index < total_R_num; ++index)
-    {
-        if (T_nonzero_num[index] != 0)
-        {
-            output_R_number++;
-        }
-    }
-
-    std::stringstream sss;
-    sss << TR_filename;
-    std::ofstream g2;
-
-    if(GlobalV::DRANK==0)
-    {
-        if (binary)
-        {
-            if(GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step)
-            {
-                g2.open(sss.str().c_str(), std::ios::binary | std::ios::app);
-            }
-            else
-            {
-                g2.open(sss.str().c_str(), std::ios::binary);
-            }
-            g2.write(reinterpret_cast<char *>(&step), sizeof(int));
-            g2.write(reinterpret_cast<char *>(&GlobalV::NLOCAL), sizeof(int));
-            g2.write(reinterpret_cast<char *>(&output_R_number), sizeof(int));
-        }
-        else
-        {
-            if(GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && step)
-            {
-                g2.open(sss.str().c_str(), std::ios::app);
-            }
-            else
-            {
-                g2.open(sss.str().c_str());
-            }
-            g2 << "STEP: " << step << std::endl;
-            g2 << "Matrix Dimension of T(R): " << GlobalV::NLOCAL <<std::endl;
-            g2 << "Matrix number of T(R): " << output_R_number << std::endl;
-        }
-    }
-
-    count = 0;
-    for (auto &R_coor : all_R_coor_ptr)
-    {
-        int dRx = R_coor.x;
-        int dRy = R_coor.y;
-        int dRz = R_coor.z;
-
-        if (T_nonzero_num[count] == 0)
-        {
-            count++;
-            continue;
-        }
-
-        if (GlobalV::DRANK == 0)
-        {
-            if (binary)
-            {
-                g2.write(reinterpret_cast<char *>(&dRx), sizeof(int));
-                g2.write(reinterpret_cast<char *>(&dRy), sizeof(int));
-                g2.write(reinterpret_cast<char *>(&dRz), sizeof(int));
-                g2.write(reinterpret_cast<char *>(&T_nonzero_num[count]), sizeof(int));
-            }
-            else
-            {
-                g2 << dRx << " " << dRy << " " << dRz << " " << T_nonzero_num[count] << std::endl;
-            }
-        }
-
-        if (GlobalV::NSPIN != 4)
-        {
-            output_single_R(g2, TR_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
-        }
-        else
-        {
-            output_soc_single_R(g2, TR_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
-        }
-
-        count++;
-
-    }
-
-    if(GlobalV::DRANK==0) 
-    {
-        g2.close();
-    }
-
-    delete[] T_nonzero_num;
-    T_nonzero_num = nullptr;
-
-    ModuleBase::timer::tick("ModuleIO","save_TR_sparse");
-    return;
-}
-
 void ModuleIO::save_dH_sparse(
     const int &istep,
     LCAO_Matrix &lm,
@@ -908,7 +622,7 @@ void ModuleIO::save_dH_sparse(
                 }
                 else
                 {
-                    output_soc_single_R(g1x[ispin], dHRx_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
+                    output_single_R(g1x[ispin], dHRx_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
                 }
             }
             if (dHy_nonzero_num[ispin][count] > 0)
@@ -919,7 +633,7 @@ void ModuleIO::save_dH_sparse(
                 }
                 else
                 {
-                    output_soc_single_R(g1y[ispin], dHRy_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
+                    output_single_R(g1y[ispin], dHRy_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
                 }
             }
             if (dHz_nonzero_num[ispin][count] > 0)
@@ -930,7 +644,7 @@ void ModuleIO::save_dH_sparse(
                 }
                 else
                 {
-                    output_soc_single_R(g1z[ispin], dHRz_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
+                    output_single_R(g1z[ispin], dHRz_soc_sparse_ptr[R_coor], sparse_threshold, binary, *lm.ParaV);
                 }
             }                        
         }
@@ -959,3 +673,119 @@ void ModuleIO::save_dH_sparse(
     ModuleBase::timer::tick("ModuleIO","save_dH_sparse");
     return;
 }
+
+template<typename Tdata>
+void ModuleIO::save_sparse(
+    const std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, Tdata>>>& smat,
+    const std::set<Abfs::Vector3_Order<int>>& all_R_coor,
+    const double& sparse_threshold,
+    const bool& binary,
+    const std::string& filename,
+    const Parallel_Orbitals& pv,
+    const std::string& label,
+    const int& istep,
+    const bool& reduce)
+{
+    ModuleBase::TITLE("ModuleIO", "save_sparse");
+    ModuleBase::timer::tick("ModuleIO", "save_sparse");
+
+    int total_R_num = all_R_coor.size();
+    std::vector<int> nonzero_num(total_R_num, 0);
+    int count = 0;
+    for (auto& R_coor : all_R_coor)
+    {
+        auto iter = smat.find(R_coor);
+        if (iter != smat.end())
+            for (auto& row_loop : iter->second)
+                nonzero_num[count] += row_loop.second.size();
+        ++count;
+    }
+    if (reduce)Parallel_Reduce::reduce_all(nonzero_num.data(), total_R_num);
+
+    int output_R_number = 0;
+    for (int index = 0; index < total_R_num; ++index)
+        if (nonzero_num[index] != 0) ++output_R_number;
+
+    std::stringstream sss;
+    sss << filename;
+    std::ofstream ofs;
+    if (!reduce || GlobalV::DRANK == 0)
+    {
+        if (binary)
+        {
+            if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && istep)
+                ofs.open(sss.str().c_str(), std::ios::binary | std::ios::app);
+            else
+                ofs.open(sss.str().c_str(), std::ios::binary);
+            ofs.write(reinterpret_cast<char*>(0), sizeof(int));
+            ofs.write(reinterpret_cast<char*>(&GlobalV::NLOCAL), sizeof(int));
+            ofs.write(reinterpret_cast<char*>(&output_R_number), sizeof(int));
+        }
+        else
+        {
+            if (GlobalV::CALCULATION == "md" && GlobalV::out_app_flag && istep)
+                ofs.open(sss.str().c_str(), std::ios::app);
+            else
+                ofs.open(sss.str().c_str());
+            ofs << "STEP: " << std::max(istep, 0) << std::endl;
+            ofs << "Matrix Dimension of " + label + "(R): " << GlobalV::NLOCAL << std::endl;
+            ofs << "Matrix number of " + label + "(R): " << output_R_number << std::endl;
+        }
+    }
+
+    count = 0;
+    for (auto& R_coor : all_R_coor)
+    {
+        int dRx = R_coor.x;
+        int dRy = R_coor.y;
+        int dRz = R_coor.z;
+
+        if (nonzero_num[count] == 0)
+        {
+            count++;
+            continue;
+        }
+
+        if (!reduce || GlobalV::DRANK == 0)
+        {
+            if (binary)
+            {
+                ofs.write(reinterpret_cast<char*>(&dRx), sizeof(int));
+                ofs.write(reinterpret_cast<char*>(&dRy), sizeof(int));
+                ofs.write(reinterpret_cast<char*>(&dRz), sizeof(int));
+                ofs.write(reinterpret_cast<char*>(&nonzero_num[count]), sizeof(int));
+            }
+            else
+            {
+                ofs << dRx << " " << dRy << " " << dRz << " " << nonzero_num[count] << std::endl;
+            }
+        }
+
+        output_single_R(ofs, smat.at(R_coor), sparse_threshold, binary, pv, reduce);
+        ++count;
+    }
+    if (!reduce || GlobalV::DRANK == 0) ofs.close();
+
+    ModuleBase::timer::tick("ModuleIO", "save_sparse");
+}
+
+template void ModuleIO::save_sparse<double>(
+    const std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>>&,
+    const std::set<Abfs::Vector3_Order<int>>&,
+    const double&,
+    const bool&,
+    const std::string&,
+    const Parallel_Orbitals&,
+    const std::string&,
+    const int&,
+    const bool&);
+template void ModuleIO::save_sparse<std::complex<double>>(
+    const std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>>&,
+    const std::set<Abfs::Vector3_Order<int>>&,
+    const double&,
+    const bool&,
+    const std::string&,
+    const Parallel_Orbitals&,
+    const std::string&,
+    const int&,
+    const bool&);
\ No newline at end of file
diff --git a/source/module_io/write_HS_sparse.h b/source/module_io/write_HS_sparse.h
index b6069c44a6..9411eeee75 100644
--- a/source/module_io/write_HS_sparse.h
+++ b/source/module_io/write_HS_sparse.h
@@ -20,25 +20,24 @@ namespace ModuleIO
         const std::string &HR_filename_up, 
         const std::string &HR_filename_down
     );
-    void save_SR_sparse(
-        LCAO_Matrix &lm,
-        const double& sparse_threshold,
-        const bool &binary,  
-        const std::string &SR_filename
-    );
-    void save_TR_sparse(
-        const int &istep,
-        LCAO_Matrix &lm,
-        const double& sparse_threshold,
-        const bool &binary,  
-        const std::string &TR_filename
-    );
     void save_dH_sparse(
         const int &istep,
         LCAO_Matrix &lm,
         const double& sparse_threshold,
-        const bool &binary
+        const bool& binary
     );
+
+    template<typename Tdata>
+    void save_sparse(
+        const std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, Tdata>>>& smat,
+        const std::set<Abfs::Vector3_Order<int>>& all_R_coor,
+        const double& sparse_threshold,
+        const bool& binary,
+        const std::string& filename,
+        const Parallel_Orbitals& pv,
+        const std::string& label,
+        const int& istep = -1,
+        const bool& reduce = true);
 }
 
 #endif
diff --git a/source/module_ri/Exx_LRI_interface.h b/source/module_ri/Exx_LRI_interface.h
index 13254909e0..9fd3e977c6 100644
--- a/source/module_ri/Exx_LRI_interface.h
+++ b/source/module_ri/Exx_LRI_interface.h
@@ -19,15 +19,22 @@ template<typename T, typename Tdata>
 class Exx_LRI_Interface
 {
 public:
+    using TC = std::array<int, 3>;
+    using TAC = std::pair<int, TC>;
+
     /// @brief  Constructor for Exx_LRI_Interface
     /// @param exx_ptr
     Exx_LRI_Interface(std::shared_ptr<Exx_LRI<Tdata>> exx_ptr) : exx_ptr(exx_ptr) {}
     Exx_LRI_Interface() = delete;
 
-    void write_Hexxs(const std::string &file_name) const;
-    void read_Hexxs(const std::string& file_name);
-    
-    using TAC = std::pair<int, std::array<int, 3>>;
+    /// read and write Hexxs using cereal
+    void write_Hexxs_cereal(const std::string& file_name) const;
+    void read_Hexxs_cereal(const std::string& file_name);
+
+    /// read and write Hexxs in CSR format
+    void write_Hexxs_csr(const std::string& file_name, const UnitCell& ucell) const;
+    void read_Hexxs_csr(const std::string& file_name, const UnitCell& ucell);
+
     std::vector< std::map<int, std::map<TAC, RI::Tensor<Tdata>>>>& get_Hexxs() const { return this->exx_ptr->Hexxs; }
     
     double& get_Eexx() const { return this->exx_ptr->Eexx; }
@@ -52,6 +59,13 @@ class Exx_LRI_Interface
         int& iter);
     int two_level_step = 0;
 private:
+
+    /// calculate CSR sparse matrix from the global matrix stored with RI::Tensor
+    /// the return type is same as LCAO_Matrix::SR_sparse,  HR_sparse, etc.
+    std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, Tdata>>>
+        calculate_RI_Tensor_sparse(const double& sparse_threshold,
+            const std::map<int, std::map<TAC, RI::Tensor<Tdata>>>& hR, const UnitCell& ucell)const;
+
     std::shared_ptr<Exx_LRI<Tdata>> exx_ptr;
     Mix_DMk_2D mix_DMk_2D;
 };
diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp
index ca13f6f977..9c04789a1a 100644
--- a/source/module_ri/Exx_LRI_interface.hpp
+++ b/source/module_ri/Exx_LRI_interface.hpp
@@ -8,28 +8,130 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h"
 
 #include <sys/time.h>
+#include "module_io/csr_reader.h"
+#include "module_io/write_HS_sparse.h"
 
 template<typename T, typename Tdata>
-void Exx_LRI_Interface<T, Tdata>::write_Hexxs(const std::string& file_name) const
+void Exx_LRI_Interface<T, Tdata>::write_Hexxs_cereal(const std::string& file_name) const
 {
-	ModuleBase::TITLE("Exx_LRI","write_Hexxs");
-	ModuleBase::timer::tick("Exx_LRI", "write_Hexxs");
-	std::ofstream ofs(file_name, std::ofstream::binary);
+    ModuleBase::TITLE("Exx_LRI", "write_Hexxs_cereal");
+    ModuleBase::timer::tick("Exx_LRI", "write_Hexxs_cereal");
+    std::ofstream ofs(file_name + "_" + std::to_string(GlobalV::MY_RANK), std::ofstream::binary);
 	cereal::BinaryOutputArchive oar(ofs);
-	oar(this->exx_ptr->Hexxs);
-	ModuleBase::timer::tick("Exx_LRI", "write_Hexxs");
+    oar(this->exx_ptr->Hexxs);
+    ModuleBase::timer::tick("Exx_LRI", "write_Hexxs_cereal");
 }
 
 template<typename T, typename Tdata>
-void Exx_LRI_Interface<T, Tdata>::read_Hexxs(const std::string& file_name)
+void Exx_LRI_Interface<T, Tdata>::read_Hexxs_cereal(const std::string& file_name)
 {
-	ModuleBase::TITLE("Exx_LRI","read_Hexxs");
-	ModuleBase::timer::tick("Exx_LRI", "read_Hexxs");
-	std::ifstream ifs(file_name, std::ofstream::binary);
+    ModuleBase::TITLE("Exx_LRI", "read_Hexxs_cereal");
+    ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal");
+    std::ifstream ifs(file_name + "_" + std::to_string(GlobalV::MY_RANK), std::ofstream::binary);
 	cereal::BinaryInputArchive iar(ifs);
 	iar(this->exx_ptr->Hexxs);
-	ModuleBase::timer::tick("Exx_LRI", "read_Hexxs");
+    ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal");
 }
+
+template<typename T, typename Tdata>
+void Exx_LRI_Interface<T, Tdata>::write_Hexxs_csr(const std::string& file_name, const UnitCell& ucell) const
+{
+    ModuleBase::TITLE("Exx_LRI", "write_Hexxs_csr");
+    ModuleBase::timer::tick("Exx_LRI", "write_Hexxs_csr");
+    std::set<Abfs::Vector3_Order<int>> all_R_coor;
+    double sparse_threshold = 1e-10;
+    for (int is = 0;is < this->exx_ptr->Hexxs.size();++is)
+    {
+        for (const auto& HexxA : this->exx_ptr->Hexxs[is])
+        {
+            const int iat0 = HexxA.first;
+            for (const auto& HexxB : HexxA.second)
+            {
+                const int iat1 = HexxB.first.first;
+                const Abfs::Vector3_Order<int> R = RI_Util::array3_to_Vector3(HexxB.first.second);
+                all_R_coor.insert(R);
+            }
+        }
+        ModuleIO::save_sparse(
+            this->calculate_RI_Tensor_sparse(sparse_threshold, this->exx_ptr->Hexxs[is], ucell),
+            all_R_coor,
+            sparse_threshold,
+            false, //binary
+            file_name + "_" + std::to_string(is) + ".csr",
+            Parallel_Orbitals(),
+            "Hexxs_" + std::to_string(is),
+            -1,
+            false);  //no reduce, one file for each process
+    }
+    ModuleBase::timer::tick("Exx_LRI", "write_Hexxs_csr");
+}
+
+template<typename T, typename Tdata>
+std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, Tdata>>>
+Exx_LRI_Interface<T, Tdata>::calculate_RI_Tensor_sparse(const double& sparse_threshold,
+    const std::map<int, std::map<TAC, RI::Tensor<Tdata>>>& hR,
+    const UnitCell& ucell)const
+{
+    ModuleBase::TITLE("Exx_LRI_Interface", "calculate_HContainer_sparse_d");
+    std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, Tdata>>> target;
+    for (auto& a1_a2R_data : hR)
+    {
+        int iat1 = a1_a2R_data.first;
+        for (auto& a2R_data : a1_a2R_data.second)
+        {
+            int iat2 = a2R_data.first.first;
+            int nw1 = ucell.atoms[ucell.iat2it[iat1]].nw;
+            int nw2 = ucell.atoms[ucell.iat2it[iat2]].nw;
+            int start1 = ucell.atoms[ucell.iat2it[iat1]].stapos_wf + ucell.iat2ia[iat1] * nw1;
+            int start2 = ucell.atoms[ucell.iat2it[iat2]].stapos_wf + ucell.iat2ia[iat2] * nw2;
+
+            const TC& R = a2R_data.first.second;
+            auto& matrix = a2R_data.second;
+            Abfs::Vector3_Order<int> dR(R[0], R[1], R[2]);
+            for (int i = 0;i < nw1;++i)
+                for (int j = 0;j < nw2;++j)
+                    target[dR][start1 + i][start2 + j] = ((std::abs(matrix(i, j)) > sparse_threshold) ? matrix(i, j) : static_cast<Tdata>(0));
+        }
+    }
+    return target;
+}
+template<typename T, typename Tdata>
+void Exx_LRI_Interface<T, Tdata>::read_Hexxs_csr(const std::string& file_name, const UnitCell& ucell)
+{
+    ModuleBase::TITLE("Exx_LRI", "read_Hexxs");
+    ModuleBase::timer::tick("Exx_LRI", "read_Hexxs");
+    this->exx_ptr->Hexxs.resize(GlobalV::NSPIN);
+    for (int is = 0;is < GlobalV::NSPIN;++is)
+    {
+        ModuleIO::csrFileReader<Tdata> csr(file_name + "_" + std::to_string(is) + ".csr");
+        int nR = csr.getNumberOfR();
+        int nbasis = csr.getMatrixDimension();
+        assert(nbasis == GlobalV::NLOCAL);
+        // allocate Hexxs[is]
+        for (int iat1 = 0; iat1 < ucell.nat; ++iat1)
+            for (int iat2 = 0;iat2 < ucell.nat;++iat2)
+                for (int iR = 0;iR < nR;++iR)
+                {
+                    const std::vector<int>& R = csr.getRCoordinate(iR);
+                    TC dR({ R[0], R[1], R[2] });
+                    this->exx_ptr->Hexxs[is][iat1][{iat2, dR}] = RI::Tensor<Tdata>({ static_cast<size_t>(ucell.atoms[ucell.iat2it[iat1]].nw), static_cast<size_t>(ucell.atoms[ucell.iat2it[iat2]].nw) });
+                }
+        // read Hexxs[is]
+        for (int i = 0;i < GlobalV::NLOCAL;++i)
+            for (int j = 0;j < GlobalV::NLOCAL;++j)
+                for (int iR = 0;iR < nR;++iR)
+                {
+                    int iat1 = ucell.iwt2iat[i];
+                    int iat2 = ucell.iwt2iat[j];
+                    const std::vector<int>& R = csr.getRCoordinate(iR);
+                    const auto& matrix = csr.getMatrix(iR);
+                    TC dR({ R[0], R[1], R[2] });
+                    this->exx_ptr->Hexxs.at(is).at(iat1).at({ iat2, dR })(ucell.iwt2iw[i], ucell.iwt2iw[j]) = matrix(i, j);
+                }
+    }
+    ModuleBase::timer::tick("Exx_LRI", "read_Hexxs");
+}
+
 template<typename T, typename Tdata>
 void Exx_LRI_Interface<T, Tdata>::exx_beforescf(const K_Vectors& kv, const Charge_Mixing& chgmix)
 {