diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp
index ba2c9c66840..8689cabef21 100644
--- a/source/module_elecstate/module_dm/density_matrix.cpp
+++ b/source/module_elecstate/module_dm/density_matrix.cpp
@@ -130,7 +130,7 @@ void DensityMatrix<TK, TR>::init_DMR(Grid_Driver* GridD_in, const UnitCell* ucel
     {
         tmp_DMR->fix_gamma();
     }
-    tmp_DMR->allocate(true);
+    tmp_DMR->allocate(nullptr, true);
     this->_DMR.push_back(tmp_DMR);
     // add another DMR if nspin==2
     if (this->_nspin == 2)
@@ -180,7 +180,7 @@ void DensityMatrix<TK, TR>::init_DMR(Record_adj& ra, const UnitCell* ucell)
     {
         tmp_DMR->fix_gamma();
     }
-    tmp_DMR->allocate(true);
+    tmp_DMR->allocate(nullptr, true);
     this->_DMR.push_back(tmp_DMR);
     // add another DMR if nspin==2
     if (this->_nspin == 2)
@@ -242,7 +242,7 @@ void DensityMatrix<TK, TR>::init_DMR(const hamilt::HContainer<TRShift>& DMR_in)
                 tmp_DMR->insert_pair(tmp_ap);
             }
         }
-        tmp_DMR->allocate(true);
+        tmp_DMR->allocate(nullptr, true);
         this->_DMR.push_back(tmp_DMR);
         if(this->_nspin == 2)
         {
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp
index a8091dc07f1..0aa7da9125a 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp
@@ -117,7 +117,7 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(Grid_Driver* Gr
     // allocate the memory of BaseMatrix in HR, and set the new values to zero
     if(std::is_same<TK, double>::value)
     {
-        this->H_V_delta->allocate(true);
+        this->H_V_delta->allocate(nullptr, true);
     }
 
     ModuleBase::timer::tick("DeePKS", "initialize_HR");
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp
index 59ca874d14c..9fea5f8eb20 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp
@@ -88,7 +88,7 @@ void hamilt::EkineticNew<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(Grid_Drive
         }
     }
     // allocate the memory of BaseMatrix in HR, and set the new values to zero
-    this->hR->allocate(true);
+    this->hR->allocate(nullptr, true);
 
     ModuleBase::timer::tick("EkineticNew", "initialize_HR");
 }
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.cpp
index fc716ae2b8d..00eb0bdf728 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.cpp
@@ -105,7 +105,7 @@ void hamilt::NonlocalNew<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(Grid_Drive
         }
     }
     // allocate the memory of BaseMatrix in HR, and set the new values to zero
-    this->hR->allocate(true);
+    this->hR->allocate(nullptr, true);
 
     ModuleBase::timer::tick("NonlocalNew", "initialize_HR");
 }
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp
index 002f44bb36b..1369bfb8dc5 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp
@@ -70,7 +70,7 @@ void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::initialize_SR(Grid_Driver
         }
     }
     // allocate the memory of BaseMatrix in SR, and set the new values to zero
-    SR->allocate(true);
+    SR->allocate(nullptr, true);
     ModuleBase::timer::tick("OverlapNew", "initialize_SR");
 }
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp
index 3a6f36530a8..a2e2b5c8598 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp
@@ -53,7 +53,7 @@ void Veff<OperatorLCAO<TK, TR>>::initialize_HR(const UnitCell* ucell_in,
         }
     }
     // allocate the memory of BaseMatrix in HR, and set the new values to zero
-    this->hR->allocate(true);
+    this->hR->allocate(nullptr, true);
 
     ModuleBase::timer::tick("Veff", "initialize_HR");
 }
diff --git a/source/module_hamilt_lcao/module_gint/gint.cpp b/source/module_hamilt_lcao/module_gint/gint.cpp
index dabdd0b693d..41edaf25b20 100644
--- a/source/module_hamilt_lcao/module_gint/gint.cpp
+++ b/source/module_hamilt_lcao/module_gint/gint.cpp
@@ -496,7 +496,7 @@ void Gint::initialize_pvpR(
 	}// end T1
 	if(npol == 1)
 	{
-		this->hRGint->allocate(0);
+		this->hRGint->allocate(nullptr, 0);
 		ModuleBase::Memory::record("Gint::hRGint",this->hRGint->get_memory_size());
 		// initialize DMRGint with hRGint when NSPIN != 4
 		for (int is = 0; is < this->DMRGint.size(); is++)
@@ -511,15 +511,15 @@ void Gint::initialize_pvpR(
 	}
 	else
 	{
-		this->hRGintCd->allocate(0);
+		this->hRGintCd->allocate(nullptr, 0);
 		ModuleBase::Memory::record("Gint::hRGintCd",this->hRGintCd->get_memory_size());
 		for (int is = 0; is < this->DMRGint.size(); is++)
 		{
-			this->DMRGint[is]->allocate(0);
+			this->DMRGint[is]->allocate(nullptr, 0);
 		}
 		ModuleBase::Memory::record("Gint::DMRGint",this->DMRGint[0]->get_memory_size() * this->DMRGint.size());
 #ifdef __MPI	
-		this->DMRGint_full->allocate(0);
+		this->DMRGint_full->allocate(nullptr, 0);
 		ModuleBase::Memory::record("Gint::DMRGint_full",this->DMRGint_full->get_memory_size());
 #endif
 	}
diff --git a/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp b/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp
index 57fc9127027..c0d3aef56ab 100644
--- a/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp
@@ -29,20 +29,17 @@ AtomPair<T>::AtomPair(const int& atom_i_, const int& atom_j_, const Parallel_Orb
     }
     this->row_size = this->paraV->get_row_size(atom_i);
     this->col_size = this->paraV->get_col_size(atom_j);
-    this->ldc = this->paraV->get_col_size();
     this->R_index.resize(3, 0);
     this->current_R = 0;
     if (existed_matrix != nullptr)
     {
-        BaseMatrix<T> tmp(row_size, col_size, (existed_matrix + row_ap * ldc + col_ap));
-        tmp.set_ldc(this->ldc);
+        BaseMatrix<T> tmp(row_size, col_size, existed_matrix);
         this->values.push_back(tmp);
     }
     else
     {
         BaseMatrix<T> tmp(row_size, col_size);
         this->values.push_back(tmp);
-        this->ldc = col_size;
     }
 }
 
@@ -65,7 +62,6 @@ AtomPair<T>::AtomPair(const int& atom_i_,
     }
     this->row_size = this->paraV->get_row_size(atom_i);
     this->col_size = this->paraV->get_col_size(atom_j);
-    this->ldc = this->paraV->get_col_size();
     this->R_index.resize(3, 0);
     this->current_R = 0;
     this->R_index[0] = rx;
@@ -73,15 +69,13 @@ AtomPair<T>::AtomPair(const int& atom_i_,
     this->R_index[2] = rz;
     if (existed_matrix != nullptr)
     {
-        BaseMatrix<T> tmp(row_size, col_size, (existed_matrix + row_ap * ldc + col_ap));
-        tmp.set_ldc(this->ldc);
+        BaseMatrix<T> tmp(row_size, col_size, existed_matrix);
         this->values.push_back(tmp);
     }
     else
     {
         BaseMatrix<T> tmp(row_size, col_size);
         this->values.push_back(tmp);
-        this->ldc = col_size;
     }
 }
 // direct save whole matrix of atom-pair
@@ -103,16 +97,13 @@ AtomPair<T>::AtomPair(const int& atom_i_,
     this->current_R = 0;
     if (existed_matrix != nullptr)
     {
-        this->ldc = row_atom_begin[natom] - row_atom_begin[0];
-        BaseMatrix<T> tmp(row_size, col_size, (existed_matrix + row_ap * ldc + col_ap));
-        tmp.set_ldc(this->ldc);
+        BaseMatrix<T> tmp(row_size, col_size, existed_matrix);
         this->values.push_back(tmp);
     }
     else
     {
         BaseMatrix<T> tmp(row_size, col_size);
         this->values.push_back(tmp);
-        this->ldc = col_size;
     }
 }
 //
@@ -140,16 +131,13 @@ AtomPair<T>::AtomPair(const int& atom_i_,
     this->R_index[2] = rz;
     if (existed_matrix != nullptr)
     {
-        this->ldc = row_atom_begin[natom] - row_atom_begin[0];
-        BaseMatrix<T> tmp(row_size, col_size, (existed_matrix + row_ap * ldc + col_ap));
-        tmp.set_ldc(this->ldc);
+        BaseMatrix<T> tmp(row_size, col_size, existed_matrix);
         this->values.push_back(tmp);
     }
     else
     {
         BaseMatrix<T> tmp(row_size, col_size);
         this->values.push_back(tmp);
-        this->ldc = col_size;
     }
 }
 
@@ -160,9 +148,8 @@ AtomPair<T>::AtomPair(const int& atom_i_, const int& atom_j_) : atom_i(atom_i_),
 
 // copy constructor
 template <typename T>
-AtomPair<T>::AtomPair(const AtomPair<T>& other)
+AtomPair<T>::AtomPair(const AtomPair<T>& other, T* data_pointer)
     : R_index(other.R_index),
-      values(other.values),
       paraV(other.paraV),
       current_R(other.current_R),
       atom_i(other.atom_i),
@@ -170,18 +157,42 @@ AtomPair<T>::AtomPair(const AtomPair<T>& other)
       row_ap(other.row_ap),
       col_ap(other.col_ap),
       row_size(other.row_size),
-      col_size(other.col_size),
-      ldc(other.ldc)
+      col_size(other.col_size)
 {
+    if(data_pointer == nullptr)
+    {
+        this->values = other.values;
+    }
+    else
+    {
+        this->values.reserve(other.values.size());
+        for(int value=0;value<other.values.size();++value)
+        {
+            hamilt::BaseMatrix<T> tmp(row_size, col_size, data_pointer);
+            this->values.push_back(tmp);
+            data_pointer += this->get_size();
+        }
+    }
 }
 
 //allocate
 template <typename T>
-void AtomPair<T>::allocate(bool is_zero)
+void AtomPair<T>::allocate(T* data_array, bool is_zero)
 {
-    for(int value=0;value<this->values.size();++value)
+    if(data_array == nullptr)
+    {
+        for(int value=0;value<this->values.size();++value)
+        {
+            this->values[value].allocate(nullptr, is_zero);
+        }
+    }
+    else
     {
-        this->values[value].allocate(is_zero);
+        for(int value=0;value<this->values.size();++value)
+        {
+            this->values[value].allocate(data_array, is_zero);
+            data_array += this->get_size();
+        }
     }
 }
 
@@ -211,7 +222,6 @@ AtomPair<T>& AtomPair<T>::operator=(const AtomPair<T>& other)
         col_ap = other.col_ap;
         row_size = other.row_size;
         col_size = other.col_size;
-        ldc = other.ldc;
     }
     return *this;
 }
@@ -228,8 +238,7 @@ AtomPair<T>::AtomPair(AtomPair<T>&& other) noexcept
       row_ap(other.row_ap),
       col_ap(other.col_ap),
       row_size(other.row_size),
-      col_size(other.col_size),
-      ldc(other.ldc)
+      col_size(other.col_size)
 {
     other.paraV = nullptr;
 }
@@ -251,7 +260,6 @@ AtomPair<T>& AtomPair<T>::operator=(AtomPair<T>&& other) noexcept
         col_ap = other.col_ap;
         row_size = other.row_size;
         col_size = other.col_size;
-        ldc = other.ldc;
     }
     return *this;
 }
@@ -367,7 +375,7 @@ BaseMatrix<T>& AtomPair<T>::get_HR_values(int rx_in, int ry_in, int rz_in)
     R_index.push_back(ry_in);
     R_index.push_back(rz_in);
     values.push_back(BaseMatrix<T>(this->row_size, this->col_size));
-    values.back().allocate(true);
+    values.back().allocate(nullptr, true);
     // return the last BaseMatrix reference in values
     return this->values.back();
 }
@@ -511,7 +519,7 @@ void AtomPair<T>::merge_to_gamma()
         {
             if(empty)
             {
-                tmp.allocate(true);
+                tmp.allocate(nullptr, true);
                 empty = false;
             }
             tmp.add_array(this->values[i].get_pointer());
@@ -520,8 +528,6 @@ void AtomPair<T>::merge_to_gamma()
     this->values.clear();
     this->values.push_back(tmp);
 
-    this->ldc = this->col_size;
-
     this->current_R = 0;
 }
 
@@ -548,7 +554,7 @@ void AtomPair<T>::add_to_matrix(std::complex<T>* hk,
             BlasConnector::axpy(this->col_size, kphase.real(), hr_tmp, 1, hk_real_pointer, 2);
             BlasConnector::axpy(this->col_size, kphase.imag(), hr_tmp, 1, hk_imag_pointer, 2);
             hk_tmp += ld_hk;
-            hr_tmp += this->ldc;
+            hr_tmp += this->col_size;
         }
     }
     // column major
@@ -562,7 +568,7 @@ void AtomPair<T>::add_to_matrix(std::complex<T>* hk,
             BlasConnector::axpy(this->col_size, kphase.real(), hr_tmp, 1, hk_real_pointer, ld_hk_2);
             BlasConnector::axpy(this->col_size, kphase.imag(), hr_tmp, 1, hk_imag_pointer, ld_hk_2);
             hk_tmp ++;
-            hr_tmp += this->ldc;
+            hr_tmp += this->col_size;
         }
     }
 }
@@ -586,23 +592,13 @@ void AtomPair<T>::add_to_matrix(T* hk, const int ld_hk, const T& kphase, const i
                 hk_tmp[nu] += matrix.get_value(mu, nu) * kphase;
             }*/
             hk_tmp += ld_hk;
-            hr_tmp += this->ldc;
+            hr_tmp += this->col_size;
         }
     }
     // column major
     else if (hk_type == 1)
     {
         hk_tmp += this->col_ap * ld_hk + this->row_ap;
-        /*for (int nu = 0; nu < this->col_size; nu++)
-        {
-            BlasConnector::axpy(this->row_size, kphase, hr_tmp, this->ldc, hk_tmp, 1);
-            for (int mu = 0; mu < this->row_size; mu++)
-            {
-                hk_tmp[mu] += matrix.get_value(mu, nu) * kphase;
-            }
-            hk_tmp += ld_hk;
-            hr_tmp++;
-        }*/
         for (int mu = 0; mu < this->row_size; mu++)
         {
             BlasConnector::axpy(this->col_size, kphase, hr_tmp, 1, hk_tmp, ld_hk);
@@ -611,7 +607,7 @@ void AtomPair<T>::add_to_matrix(T* hk, const int ld_hk, const T& kphase, const i
                 hk_tmp[mu] += matrix.get_value(mu, nu) * kphase;
             }*/
             ++hk_tmp;
-            hr_tmp += this->ldc;
+            hr_tmp += this->col_size;
         }
     }
 }
@@ -639,18 +635,10 @@ void AtomPair<T>::add_to_array(std::complex<T>* array, const std::complex<T>& kp
 }
 
 template <typename T>
-T& AtomPair<T>::get_matrix_value(const size_t& i_row_global, const size_t& j_col_global) const
+std::tuple<std::vector<int>, T*> AtomPair<T>::get_matrix_values(int ir) const
 {
-    int i_row_local = this->paraV == nullptr ? i_row_global : this->paraV->global2local_row(i_row_global);
-    int j_col_local = this->paraV == nullptr ? j_col_global : this->paraV->global2local_col(j_col_global);
-#ifdef __DEBUG
-    assert(i_row_local != -1 && j_col_local != -1);
-    assert(current_R < this->values.size());
-    assert(current_R >= 0);
-#endif
-    size_t i_row_in = i_row_local - row_ap;
-    size_t j_col_in = j_col_local - col_ap;
-    return this->values[current_R].get_value(i_row_in, j_col_in);
+    if(ir<0) ir = this->current_R;
+    return std::tuple<std::vector<int>, T*>({this->row_ap, this->row_size, this->col_ap, this->col_size}, this->values[ir].get_pointer());
 }
 
 // interface for get (rx, ry, rz) of index-th R-index in this->R_index, the return should be int[3]
diff --git a/source/module_hamilt_lcao/module_hcontainer/atom_pair.h b/source/module_hamilt_lcao/module_hcontainer/atom_pair.h
index d13650f6e7a..86539891d12 100644
--- a/source/module_hamilt_lcao/module_hcontainer/atom_pair.h
+++ b/source/module_hamilt_lcao/module_hcontainer/atom_pair.h
@@ -82,7 +82,7 @@ class AtomPair
              T* existed_matrix = nullptr);
 
     // copy constructor
-    AtomPair(const AtomPair<T>& other);
+    AtomPair(const AtomPair<T>& other, T* data_pointer = nullptr);
     // move constructor
     AtomPair(AtomPair&& other) noexcept;
 
@@ -96,7 +96,7 @@ class AtomPair
     /**
      * @brief allocate memory for all the BaseMatrix
     */
-    void allocate(bool if_zero = false);
+    void allocate(T* data_array = nullptr, bool if_zero = false);
 
     /**
      * @brief set values in every BaseMatrix to zero
@@ -176,7 +176,19 @@ class AtomPair
     // these four interface can be used only when R-index has been choosed (current_R >= 0)
     T& get_value(const int& i) const;
     T& get_value(const int& row, const int& col) const;
-    T& get_matrix_value(const size_t& i_row_global, const size_t& j_col_global) const;
+
+    /**
+     * @brief get values of this->values[ir] for a whole matrix
+     * @param ir index of this->values
+     * @return std::tuple<std::vector<int>, T*>
+     * std::vector<int>(4) contains (row_begin_index, row_size, col_begin_index, col_size)
+     * T* is pointer of values[ir].value_begin, legal index is [0, row_size*col_size)
+    */
+    std::tuple<std::vector<int>, T*> get_matrix_values(int ir = -1) const;
+
+    /**
+     * @brief get pointer of value from a submatrix
+    */
     T* get_pointer(int ir=-1) const;
 
     // add another BaseMatrix<T> to this->values with specific R index.
@@ -270,7 +282,6 @@ class AtomPair
     int col_ap = -1;
     int row_size = 0;
     int col_size = 0;
-    int ldc = -1; // leading_dimention_colomn
 };
 
 } // namespace hamilt
diff --git a/source/module_hamilt_lcao/module_hcontainer/base_matrix.cpp b/source/module_hamilt_lcao/module_hcontainer/base_matrix.cpp
index 1a7b80d7b50..ff02d7744ab 100644
--- a/source/module_hamilt_lcao/module_hcontainer/base_matrix.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/base_matrix.cpp
@@ -12,20 +12,8 @@ BaseMatrix<T>::BaseMatrix(const int& nrow_, const int& ncol_, T* data_existed)
 {
     nrow_local = nrow_;
     ncol_local = ncol_;
-    // the default memory_type is 1 (dense matrix), it doesn't matter for initialization
-    if (data_existed == nullptr)
-    {
-        this->allocated = false;
-        this->ldc = ncol_local;
-        this->memory_type = 1;
-    }
-    else
-    {
-        value_begin = data_existed;
-        this->allocated = false;
-        this->ldc = ncol_local;
-        this->memory_type = 2;
-    }
+    this->allocated = false;
+    value_begin = data_existed;
 }
 
 // move constructor
@@ -49,8 +37,6 @@ BaseMatrix<T>::BaseMatrix(const BaseMatrix<T>& matrix)
 {
     this->nrow_local = matrix.nrow_local;
     this->ncol_local = matrix.ncol_local;
-    this->memory_type = matrix.memory_type;
-    this->ldc = matrix.ldc;
     if (matrix.allocated)
     {
         this->value_begin = new T[nrow_local * ncol_local];
@@ -79,52 +65,44 @@ BaseMatrix<T>::~BaseMatrix()
 
 // allocate
 template <typename T>
-void BaseMatrix<T>::allocate(bool if_zero)
+void BaseMatrix<T>::allocate(T* data_array, bool if_zero)
 {
 #ifdef __DEBUG
 assert(nrow_local*ncol_local>0);
 #endif
-    if(this->value_begin == nullptr)
+    if(data_array != nullptr && !this->allocated)
     {
-        this->value_begin = new T[nrow_local * ncol_local];
-        if(if_zero) 
-        {
-            this->set_zero();
-        }
-        this->allocated = true;
+        this->value_begin = data_array;
     }
-}
-
-// zeros
-template <typename T>
-void BaseMatrix<T>::set_zero()
-{
-#ifdef __DEBUG
-assert(this->value_begin != nullptr);
-#endif
-    if(this->memory_type == 1)
+    else if(data_array != nullptr && this->allocated)
     {
-        ModuleBase::GlobalFunc::ZEROS(this->value_begin, nrow_local * ncol_local);
+        delete[] this->value_begin;
+        this->value_begin = data_array;
+        this->allocated = false;
     }
-    else if(this->memory_type == 2)
+    else if(data_array == nullptr && !this->allocated)
     {
-        for(int i = 0; i < nrow_local; i++)
-        {
-            ModuleBase::GlobalFunc::ZEROS(this->value_begin + i * this->ldc, ncol_local);
-        }
+        this->value_begin = new T[nrow_local * ncol_local];
+        this->allocated = true;
     }
     else
     {
-        std::cout << "Error: memory_type is not defined!" << std::endl;
-        exit(1);
+        // do nothing
+    }
+    if(if_zero) 
+    {
+        this->set_zero();
     }
 }
 
-// set_memory_type
+// zeros
 template <typename T>
-void BaseMatrix<T>::set_ldc(const int& ldc_in)
+void BaseMatrix<T>::set_zero()
 {
-    this->ldc = ldc_in;
+#ifdef __DEBUG
+assert(this->value_begin != nullptr);
+#endif
+    ModuleBase::GlobalFunc::ZEROS(this->value_begin, nrow_local * ncol_local);
 }
 
 // add_array
@@ -133,42 +111,23 @@ void BaseMatrix<T>::add_array(T* array)
 {
     // if allocated, save data from array into matrix
     // if whole matrix and 2d-block format, save data from array into matrix either
-    if (this->allocated || this->memory_type == 0 || this->memory_type == 1)
+    for (int i = 0; i < nrow_local * ncol_local; ++i)
     {
-        for (int i = 0; i < nrow_local * ncol_local; i++)
-        {
-            value_begin[i] += array[i];
-        }
-    }
-    else
-    { // if not allocated, then it is a wrapper of block submatrix
-        if (this->memory_type == 2)
-        {
-            for (int i = 0; i < this->nrow_local; i++)
-            {
-                for (int j = 0; j < this->ncol_local; j++)
-                {
-                    this->add_element(i, j, *array);
-                    array++;
-                }
-            }
-        }
+        value_begin[i] += array[i];
     }
 }
 
 template <typename T>
 void BaseMatrix<T>::add_element(int mu, int nu, const T& value)
 {
-    const int dim = this->memory_type == 2 ? this->ldc : this->ncol_local;
-    int index = mu * dim + nu;
+    int index = mu * this->ncol_local + nu;
     value_begin[index] += value;
 }
 
 template <typename T>
 T& BaseMatrix<T>::get_value(const size_t& i_row, const size_t& j_col) const
 {
-    const int dim = this->memory_type == 2 ? this->ldc : this->ncol_local;
-    int index = i_row * dim + j_col;
+    int index = i_row * this->ncol_local + j_col;
     return value_begin[index];
 }
 
@@ -186,14 +145,12 @@ BaseMatrix<T>& BaseMatrix<T>::operator=(const BaseMatrix<T>& other)
     {
         this->nrow_local = other.nrow_local;
         this->ncol_local = other.ncol_local;
-        this->memory_type = other.memory_type;
-        this->ldc = other.ldc;
         if (other.allocated)
         {
             this->value_begin = new T[nrow_local * ncol_local];
             ModuleBase::GlobalFunc::ZEROS(this->value_begin, nrow_local * ncol_local);
             this->allocated = true;
-            for (int i = 0; i < nrow_local * ncol_local; i++)
+            for (int i = 0; i < nrow_local * ncol_local; ++i)
             {
                 this->value_begin[i] = other.value_begin[i];
             }
diff --git a/source/module_hamilt_lcao/module_hcontainer/base_matrix.h b/source/module_hamilt_lcao/module_hcontainer/base_matrix.h
index de55ccc30c9..d866c6f67fe 100644
--- a/source/module_hamilt_lcao/module_hcontainer/base_matrix.h
+++ b/source/module_hamilt_lcao/module_hcontainer/base_matrix.h
@@ -12,11 +12,6 @@ namespace hamilt
  * It has two ways to arrange data:
  * 1. allocate data itself
  * 2. only access data but be arranged by other class
- * It has four ways to access data:
- * 1. whole matrix
- * 2. 2d-block
- * 3. submatrix in whole matrix
- * 4. sparse matrix in whole matrix , not implemented yet
  */
 template <typename T>
 class BaseMatrix
@@ -36,12 +31,10 @@ class BaseMatrix
      * if this->value_begin is not nullptr, it will be neglected
      * if this->value_begin is nullptr, it will allocate memory with size nrow_local * ncol_local
     */
-    void allocate(bool if_zero = false);
+    void allocate(T* data_array = nullptr, bool if_zero = false);
 
     /**
      * @brief set value in the matrix to zero
-     * if memory_type == 1 , it will set all value in the matrix to zero
-     * if memory_type == 2 , it will set all value in the submatrix to zero
     */
     void set_zero();
 
@@ -62,8 +55,6 @@ class BaseMatrix
     // for inside matrix
     /**
      * @brief get value from a whole matrix
-     * for memory_type = 0 or 1, ncol_local will be used to calculate the index
-     * for memory_type = 2, ldc will be used to calculate the index
      *
      * @param i_row row index
      * @param j_col column index
@@ -75,8 +66,6 @@ class BaseMatrix
      */
     T* get_pointer() const;
 
-    void set_ldc(const int& ldc_in);
-
     // operator= for copy assignment
     BaseMatrix& operator=(const BaseMatrix& other);
 
@@ -102,15 +91,6 @@ class BaseMatrix
     // number of rows and columns
     int nrow_local = 0;
     int ncol_local = 0;
-
-    // memory type, choose how to access value via pointer of array
-    // 1 is dense matrix
-    // 2 is submatrix in whole matrix
-    // 3 is sparse matrix in whole matrix , not implemented yet
-    int memory_type = 1;
-
-    // leading dimension of matrix, used with memory_type = 2
-    int ldc = 0;
 };
 
 } // namespace hamilt
diff --git a/source/module_hamilt_lcao/module_hcontainer/func_transfer.cpp b/source/module_hamilt_lcao/module_hcontainer/func_transfer.cpp
index 312857cf274..1da0c524a18 100644
--- a/source/module_hamilt_lcao/module_hcontainer/func_transfer.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/func_transfer.cpp
@@ -586,6 +586,49 @@ void transferParallels2Serials(const hamilt::HContainer<TR>& hR_p, hamilt::HCont
 
 }
 
+template<typename TR>
+void gatherParallels(const hamilt::HContainer<TR>& hR_p,
+                     hamilt::HContainer<TR>* hR_s,
+                     const int serial_rank)
+{
+    // gather <IJR>s from all ranks to serial_rank
+    int my_rank, size;
+    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &size);
+    std::vector<int> para_ijrs = hR_p.get_ijr_info();
+    if (my_rank == serial_rank)
+    {
+        hR_s->insert_ijrs(&para_ijrs);
+        for (int i = 0; i < size; ++i)
+        {
+            if (i == serial_rank)
+                continue;
+            std::vector<int> tmp_ijrs;
+            MPI_Status status;
+            int tmp_size = 0;
+            MPI_Recv(&tmp_size, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &status);
+            tmp_ijrs.resize(tmp_size);
+            MPI_Recv(tmp_ijrs.data(),
+                     tmp_ijrs.size(),
+                     MPI_INT,
+                     i,
+                     0,
+                     MPI_COMM_WORLD,
+                     &status);
+            hR_s->insert_ijrs(&tmp_ijrs);
+        }
+        hR_s->allocate();
+    }
+    else
+    {
+        int tmp_size = para_ijrs.size();
+        MPI_Send(&tmp_size, 1, MPI_INT, serial_rank, 0, MPI_COMM_WORLD);
+        MPI_Send(para_ijrs.data(), para_ijrs.size(), MPI_INT, serial_rank, 0, MPI_COMM_WORLD);
+    }
+    // gather values from Parallels to target serial_rank
+    transferParallels2Serial(hR_p, hR_s, serial_rank);
+}
+
 // specialize for double and std::complex<double>
 template void transferSerial2Parallels(const hamilt::HContainer<double>& hR_s,
                                       hamilt::HContainer<double>* hR_p,
@@ -608,6 +651,12 @@ template void transferParallels2Serials(const hamilt::HContainer<double>& hR_p,
                                         hamilt::HContainer<double>* hR_s);
 template void transferParallels2Serials(const hamilt::HContainer<std::complex<double>>& hR_p,
                                         hamilt::HContainer<std::complex<double>>* hR_s);
+template void gatherParallels(const hamilt::HContainer<double>& hR_p,
+                              hamilt::HContainer<double>* hR_s,
+                              const int serial_rank);
+template void gatherParallels(const hamilt::HContainer<std::complex<double>>& hR_p,
+                                hamilt::HContainer<std::complex<double>>* hR_s,
+                                const int serial_rank);
 
 } // namespace hamilt
 
diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp b/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp
index 1984d1049f5..c81e9e92873 100644
--- a/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp
@@ -15,14 +15,29 @@ HContainer<T>::~HContainer()
 
 // copy constructor
 template <typename T>
-HContainer<T>::HContainer(const HContainer<T>& HR_in)
+HContainer<T>::HContainer(const HContainer<T>& HR_in, T* data_array)
 {
-    this->atom_pairs = HR_in.atom_pairs;
     this->sparse_ap = HR_in.sparse_ap;
     this->sparse_ap_index = HR_in.sparse_ap_index;
     this->gamma_only = HR_in.gamma_only;
     this->paraV = HR_in.paraV;
     this->current_R = -1;
+    this->wrapper_pointer = data_array;
+    if(data_array == nullptr)
+    {
+        this->atom_pairs = HR_in.atom_pairs;
+    }
+    else
+    {
+        this->atom_pairs.reserve(HR_in.atom_pairs.size());
+        for(int iap=0;iap<HR_in.atom_pairs.size();iap++)
+        {
+            hamilt::AtomPair<T>& target = HR_in.get_atom_pair(iap);
+            hamilt::AtomPair<T> tmp(target, data_array);
+            data_array += target.get_R_size() * target.get_size();
+            this->atom_pairs.push_back(tmp);
+        }
+    }
     // tmp terms not copied
 }
 
@@ -107,42 +122,58 @@ HContainer<T>::HContainer(const UnitCell& ucell_, const Parallel_Orbitals* paraV
             }
         }
     }
-    this->allocate(true);
+    this->allocate(nullptr, true);
 }
 
 //HContainer(const Parallel_Orbitals* paraV, T* data_pointer = nullptr);
 template <typename T>
-HContainer<T>::HContainer(const Parallel_Orbitals* paraV_in, T* data_pointer)
+HContainer<T>::HContainer(const Parallel_Orbitals* paraV_in, T* data_pointer, const std::vector<int>* ijr_info)
 {
     this->current_R = -1;
-    // use HContainer as a wrapper
-    if(data_pointer != nullptr)
-    {
-        this->gamma_only = true;
-        this->wrapper_pointer = data_pointer;
-    }
-    else // use HContainer as a container
-    {
-        this->gamma_only = false;
-    }
+
+    // use HContainer as a wrapper(!nullptr) or container(nullptr)
+    this->wrapper_pointer = data_pointer;
+
+#ifdef __DEBUG
+    assert(paraV_in != nullptr);
+#endif
+
     // save Parallel_Orbitals pointer
     this->paraV = paraV_in;
     // initialize sparse_ap
     int natom = paraV->atom_begin_row.size();
     this->sparse_ap.resize(natom);
     this->sparse_ap_index.resize(natom);
+    if(ijr_info != nullptr)
+    {
+        this->insert_ijrs(ijr_info);
+        // allocate memory
+        this->allocate(data_pointer, false);
+    }
 }
 
 // allocate
 template <typename T>
-void HContainer<T>::allocate(bool is_zero)
+void HContainer<T>::allocate(T* data_array, bool is_zero)
 {
+    if(data_array == nullptr)
+    {
 #ifdef _OPENMP
 #pragma omp parallel for
 #endif
-    for(int it=0;it<this->atom_pairs.size();it++)
+        for(int it=0;it<this->atom_pairs.size();it++)
+        {
+            this->atom_pairs[it].allocate(nullptr, is_zero);
+        }
+    }
+    else
     {
-        this->atom_pairs[it].allocate(is_zero);
+        for(int it=0;it<this->atom_pairs.size();it++)
+        {
+            this->atom_pairs[it].allocate(data_array, is_zero);
+            // move data_array pointer for the next AtomPair
+            data_array += this->atom_pairs[it].get_R_size() * this->atom_pairs[it].get_size();
+        }
     }
 }
 
@@ -568,6 +599,25 @@ size_t HContainer<T>::get_memory_size() const
     return memory;
 }
 
+// get_nnr
+template <typename T>
+size_t HContainer<T>::get_nnr() const
+{
+    size_t sum = 0;
+    for(int iap=0;iap < this->atom_pairs.size();++iap)
+    {
+        sum += this->atom_pairs[iap].get_R_size() * this->atom_pairs[iap].get_size();
+    }
+    return sum;
+}
+
+// get_wrapper
+template <typename T>
+T* HContainer<T>::get_wrapper() const
+{
+    return this->wrapper_pointer;
+}
+
 // synchronize
 template <typename T>
 void HContainer<T>::shape_synchron( const HContainer<T>& other)
@@ -608,6 +658,69 @@ void HContainer<T>::shape_synchron( const HContainer<T>& other)
     }
 }
 
+// get_IJR_info
+template <typename T>
+std::vector<int> HContainer<T>::get_ijr_info() const
+{
+    // get number of atom pairs
+    std::vector<int> ijr_info;
+    ijr_info.push_back(this->atom_pairs.size());
+    // loop atom pairs
+    for (int i = 0; i < this->atom_pairs.size(); ++i)
+    {
+        // get atom_i and atom_j
+        const int atom_i = this->atom_pairs[i].get_atom_i();
+        const int atom_j = this->atom_pairs[i].get_atom_j();
+        // push back atom_i, atom_j
+        ijr_info.push_back(atom_i);
+        ijr_info.push_back(atom_j);
+        // get number of R
+        const int number_R = this->atom_pairs[i].get_R_size();
+        ijr_info.push_back(number_R);
+        // loop R
+        for (int ir = 0; ir < number_R; ++ir)
+        {
+            int* R_pointer = this->atom_pairs[i].get_R_index(ir);
+            ijr_info.push_back(R_pointer[0]);
+            ijr_info.push_back(R_pointer[1]);
+            ijr_info.push_back(R_pointer[2]);
+        }
+    }
+    return ijr_info;
+}
+
+template<typename T>
+void HContainer<T>::insert_ijrs(const std::vector<int>* ijrs)
+{
+    if (this->paraV == nullptr)
+    {
+        ModuleBase::WARNING_QUIT("HContainer::insert_ijrs", "paraV pointer can not be nullptr!");
+    }
+    // get number of atom pairs
+    const int number_ap = (*ijrs)[0];
+    // loop AtomPairs and unpack values
+    const int* ijr_p = ijrs->data() + 1;
+    for (int i = 0; i < number_ap; ++i)
+    {
+        // get atom_i and atom_j
+        const int atom_i = *ijr_p++;
+        const int atom_j = *ijr_p++;
+        // get number of R
+        const int number_R = *ijr_p++;
+        for (int k = 0; k < number_R; ++k)
+        {
+            int r_index[3];
+            // get R index
+            r_index[0] = *ijr_p++;
+            r_index[1] = *ijr_p++;
+            r_index[2] = *ijr_p++;
+            //insert this IJ AtomPair
+            AtomPair<T> tmp_ap(atom_i, atom_j, r_index[0], r_index[1], r_index[2], this->paraV);
+            this->insert_pair(tmp_ap);
+        }
+    }
+}
+
 // T of HContainer can be double or complex<double>
 template class HContainer<double>;
 template class HContainer<std::complex<double>>;
diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer.h b/source/module_hamilt_lcao/module_hcontainer/hcontainer.h
index 567a9ce8e36..afbf636ed85 100644
--- a/source/module_hamilt_lcao/module_hcontainer/hcontainer.h
+++ b/source/module_hamilt_lcao/module_hcontainer/hcontainer.h
@@ -165,8 +165,11 @@ class HContainer
     // Destructor of class HContainer
     ~HContainer();
 
-    // copy constructor
-    HContainer(const HContainer<T>& HR_in);
+    /**
+     * @brief copy constructor
+     * when data_array is not nullptr, new HContainer will be wrapper for data_array
+    */
+    HContainer(const HContainer<T>& HR_in, T* data_array = nullptr);
 
     // move constructor
     HContainer(HContainer<T>&& HR_in);
@@ -183,13 +186,21 @@ class HContainer
      * it will not allocate memory for atom_pairs
      * this case will forbit inserting empty atom_pair
     */
-    HContainer(const Parallel_Orbitals* paraV, T* data_pointer = nullptr);
+    HContainer(
+      const Parallel_Orbitals* paraV, 
+      T* data_pointer = nullptr, 
+      const std::vector<int>* ijr_info = nullptr);
 
     /**
      * @brief allocate memory for all <IJR> matrix
+     * if data_array is not nullptr, 
+     *     use memory after data_array for each BaseMatrix;
+     *     if BaseMatrix has memory allocated before, it will be freed first. 
+     * if data_array is nullptr, allocate memory for each BaseMatrix
+     * @param data_array pointer of data array
      * @param if_zero if true, set all values to zero
     */
-    void allocate(bool if_zero = false);
+    void allocate(T* data_array = nullptr, bool if_zero = false);
 
     /**
      * @brief set values of all <IJR> matrix to zero
@@ -357,6 +368,32 @@ class HContainer
     */
     size_t get_memory_size() const;
 
+    /**
+     * @brief calculate total size of data in HContainer, 
+     * named nnr inherited from history
+     * all AtomPairs and BaseMatrixs are counted
+    */
+    size_t get_nnr() const;
+
+    /**
+     * @brief get infomation of IJR pairs in HContainer
+     * the return vector format is {size_I, I1, size_J, J1, size_R, R1x, R1y, R1z, ..., J2, ...}
+    */
+    std::vector<int> get_ijr_info() const;
+
+    /**
+     * @brief use infomation of IJ pairs to expand HContainer
+     * the input vector format is {size_IJ_pairs, I1, J1, size_R, R1x, R1y, R1z, ..., I2, J2, ...}
+     * HContainer has not been allocated after this function, 
+     * user should call allocate(...) to allocate memory.
+    */
+    void insert_ijrs(const std::vector<int>* ijrs);
+
+    /**
+     * @brief return the wrapper_pointer
+    */
+    T* get_wrapper() const;
+
     /**
      * @brief synchronization of atom-pairs for read-in HContainer
      * new <IJR> pair from read-in HContainer will be inserted into this->atom-pairs
diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h b/source/module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h
index 644b7d071a0..31bf65e3da0 100644
--- a/source/module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h
+++ b/source/module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h
@@ -65,6 +65,18 @@ template<typename TR>
 void transferParallels2Serials(const hamilt::HContainer<TR>& hR_p,
                              hamilt::HContainer<TR>* hR_s);
 
+/**
+ * @brief gather the HContainer from all parallel objects to target serial object
+ * the serial object should be empty before gather
+ * @param hR_p the HContainer of <I,J,R> atom pairs in parallel object
+ * @param hR_s the empty HContainer of <I,J,R> atom pairs in serial object
+ * @param serial_rank the rank of target serial object
+*/
+template<typename TR>
+void gatherParallels(const hamilt::HContainer<TR>& hR_p,
+                     hamilt::HContainer<TR>* hR_s,
+                     const int serial_rank);
+
 #endif // __MPI
 
 } // namespace hamilt
\ No newline at end of file
diff --git a/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp b/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp
index 9a2bbd570c3..c2df1b4d777 100644
--- a/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp
@@ -70,24 +70,15 @@ void Output_HContainer<T>::write_single_R(int rx, int ry, int rz)
     for (int iap = 0; iap < this->_hcontainer->size_atom_pairs(); ++iap)
     {
         auto atom_pair = this->_hcontainer->get_atom_pair(iap);
-        int iat1 = atom_pair.get_atom_i();
-        int iat2 = atom_pair.get_atom_j();
-        int T1 = _ucell.iat2it[iat1];
-        int T2 = _ucell.iat2it[iat2];
-        int I1 = _ucell.iat2ia[iat1];
-        int I2 = _ucell.iat2ia[iat2];
-        const int start1 = _ucell.itiaiw2iwt(T1, I1, 0);
-        const int start2 = _ucell.itiaiw2iwt(T2, I2, 0);
-        int size1 = _ucell.atoms[T1].nw;
-        int size2 = _ucell.atoms[T2].nw;
-        for (int iw1 = 0; iw1 < size1; ++iw1)
+        auto tmp_matrix_info = atom_pair.get_matrix_values();
+        int* tmp_index = std::get<0>(tmp_matrix_info).data();
+        T* tmp_data = std::get<1>(tmp_matrix_info);
+        for (int irow = tmp_index[0]; irow < tmp_index[0] + tmp_index[1]; ++irow)
         {
-            const int global_index1 = start1 + iw1;
-            for (int iw2 = 0; iw2 < size2; ++iw2)
+            for (int icol = tmp_index[2]; icol < tmp_index[2] + tmp_index[3]; ++icol)
             {
-                const int global_index2 = start2 + iw2;
-                T tmp_matrix_value = atom_pair.get_matrix_value(global_index1, global_index2);
-                sparse_matrix.insert(global_index1, global_index2, tmp_matrix_value);
+                sparse_matrix.insert(irow, icol, *tmp_data);
+                tmp_data++;
                 // to do: consider 2D block-cyclic distribution
             }
         }
diff --git a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer.cpp b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer.cpp
index 0a818256497..26241825567 100644
--- a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer.cpp
@@ -370,10 +370,9 @@ TEST_F(HContainerTest, data)
 
 // using TEST_F to test functions in BaseMatrix
 // 1. test constructor with existed data
-// 2. test set_ldc
-// 3. test add_array with memory_type = 2
-// 4. test add_element
-// 5. test get_value
+// 2. test add_array with memory_type = 2
+// 3. test add_element
+// 4. test get_value
 TEST_F(HContainerTest, basematrix_funcs)
 {
     // 1. test constructor with existed data
@@ -393,38 +392,7 @@ TEST_F(HContainerTest, basematrix_funcs)
     EXPECT_EQ(BM_copy.get_value(0, 1), 2);
     EXPECT_EQ(BM_copy.get_value(1, 0), 3);
     EXPECT_EQ(BM_copy.get_value(1, 1), 4);
-    // 2. test set_ldc
-    BM.set_ldc(3);
-    EXPECT_EQ(BM.get_value(0, 0), 1);
-    EXPECT_EQ(BM.get_value(1, 0), 4);
-    BM.set_ldc(2);
-    EXPECT_EQ(BM.get_value(1, 0), 3);
-    // 3. test add_array with memory_type = 2
-    double data_ptr2[4] = {5, 6, 7, 8};
-    BM.add_array(&data_ptr2[0]);
-    // check if data is correct
-    EXPECT_EQ(BM.get_value(0, 0), 6);
-    EXPECT_EQ(data_ptr[0], 6);
-    EXPECT_EQ(BM.get_value(0, 1), 8);
-    EXPECT_EQ(BM.get_value(1, 0), 10);
-    EXPECT_EQ(BM.get_value(1, 1), 12);
-    // 4. test add_element
-    BM.add_element(0, 0, 1);
-    BM.add_element(0, 1, 2);
-    BM.add_element(1, 0, 3);
-    BM.add_element(1, 1, 4);
-    // check if data is correct
-    EXPECT_EQ(BM.get_value(0, 0), 7);
-    EXPECT_EQ(BM.get_value(0, 1), 10);
-    EXPECT_EQ(BM.get_value(1, 0), 13);
-    EXPECT_EQ(BM.get_value(1, 1), 16);
-    // 5. test get_value
-    hamilt::BaseMatrix<double> BM2(2, 2);
-    BM2.allocate(true);
-    BM2.add_element(0, 0, 1);
-    BM2.add_element(0, 1, 2);
-    EXPECT_EQ(BM2.get_value(0, 0), 1);
-    EXPECT_EQ(BM2.get_value(0, 1), 2);
+    
 } 
 
 // using TEST_F to test functions in AtomPair
@@ -490,7 +458,7 @@ TEST_F(HContainerTest, atompair_funcs)
         for(int atom_j = 0; atom_j<2; ++atom_j)
         {
             hamilt::AtomPair<double> tmp(atom_i, atom_j, 0, 0, 0, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, nullptr);
-            tmp.allocate(false);
+            tmp.allocate(nullptr, false);
             double* tmp_data = tmp.get_HR_values(0, 0, 0).get_pointer();
             for(int i=0;i<4;++i)
             {
@@ -554,7 +522,6 @@ TEST_F(HContainerTest, atompair_funcs)
             EXPECT_EQ(hk_data3[j*4+i], hk_data2_correct[i*4+j]);
         }
     }
-
     // 6. add_to_array
     std::vector<std::complex<double>> hr_array(16, 0.0);
     std::vector<double> hr_array2(16, 0.0);
@@ -578,7 +545,12 @@ TEST_F(HContainerTest, atompair_funcs)
     // check hr_array and hr_array2 are correct
     std::complex<double> correct1;
     double correct2;
-    double correct_array[16] = {1, 2, 3, 4, 5, 6, 7, 8, 3, 4, 5, 6, 7, 8, 9, 10};
+    double correct_array[16] = {
+        1, 2, 3, 4, 
+        5, 6, 7, 8, 
+        3, 4, 5, 6, 
+        7, 8, 9, 10};
+    double test_array[16] = {1, 2, 5, 6, 3, 4, 7, 8, 3, 4, 7, 8, 5, 6, 9, 10};
     for(int i=0;i<4;++i)
     {
         for(int j=0;j<4;++j)
@@ -590,33 +562,39 @@ TEST_F(HContainerTest, atompair_funcs)
         }
     }
     // construct AtomPair from existed matrix
-    hamilt::AtomPair<double> atom_ij4(0, 0, &PO, correct_array);
-    EXPECT_EQ(atom_ij4.get_value(0, 0), correct_array[0]);
-    EXPECT_EQ(atom_ij4.get_value(1, 1), correct_array[5]);
-    EXPECT_EQ(atom_ij4.get_value(0), correct_array[0]);
-    hamilt::AtomPair<double> atom_ij5(0, 1, 1, 1, 1, &PO, correct_array);
-    hamilt::AtomPair<double> atom_ij6(1, 0, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, correct_array);
-    hamilt::AtomPair<double> atom_ij7(1, 1, 1, 1, 1, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, correct_array);
+    hamilt::AtomPair<double> atom_ij4(0, 0, &PO, test_array);
+    EXPECT_EQ(atom_ij4.get_value(0, 0), test_array[0]);
+    EXPECT_EQ(atom_ij4.get_value(1, 1), test_array[3]);
+    EXPECT_EQ(atom_ij4.get_value(0), test_array[0]);
+    hamilt::AtomPair<double> atom_ij5(0, 1, 1, 1, 1, &PO, &test_array[4]);
+    hamilt::AtomPair<double> atom_ij6(1, 0, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, &test_array[8]);
+    hamilt::AtomPair<double> atom_ij7(1, 1, 1, 1, 1, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, &test_array[12]);
     // get_matrix_value will use global2local_row and global2local_col in Parallel_Orbitals
     // so we need to set them
     std::ofstream ofs("test_hcontainer.log");
     PO.set_global2local(4, 4, false, ofs);
-    for(int i=0;i<4;++i)
-    {
-        for(int j=0;j<4;++j)
+    // joint 4 2*2 BaseMatrix to whole 4*4 matrix
+    // lambda function for check data
+    auto checkdata = [&](hamilt::AtomPair<double>& ap_in) {
+        auto data_ij4 = ap_in.get_matrix_values();
+        int* tmp_index = std::get<0>(data_ij4).data();
+        double* tmp_data = std::get<1>(data_ij4);
+        double sum_error = 0.0;
+        for(int irow = tmp_index[0]; irow < tmp_index[0] + tmp_index[1]; ++irow)
         {
-            correct2 = correct_array[i*4+j];
-            hamilt::AtomPair<double>* atom_ij_ptr; 
-            if(i/2==0 && j/2==0){ atom_ij_ptr = &atom_ij4; }
-            else if(i/2==0 && j/2==1){ atom_ij_ptr = &atom_ij5; }
-            else if(i/2==1 && j/2==0){ atom_ij_ptr = &atom_ij6; }
-            else if(i/2==1 && j/2==1){ atom_ij_ptr = &atom_ij7; }
-            else{ atom_ij_ptr = nullptr; }
-            EXPECT_EQ(atom_ij_ptr->get_matrix_value(i, j), correct2);
+            for(int icol = tmp_index[2]; icol < tmp_index[2] + tmp_index[3]; ++icol)
+            {
+                sum_error += std::abs(*tmp_data++ - correct_array[irow*4+icol]);
+            }
         }
-    }
+        return sum_error;
+    };
+    EXPECT_EQ(checkdata(atom_ij4), 0.0);
+    EXPECT_EQ(checkdata(atom_ij5), 0.0);
+    EXPECT_EQ(checkdata(atom_ij6), 0.0);
+    EXPECT_EQ(checkdata(atom_ij7), 0.0);
 
-    hamilt::HContainer<double> HR_wrapper(&PO, correct_array);
+    hamilt::HContainer<double> HR_wrapper(&PO, test_array);
     HR_wrapper.insert_pair(atom_ij4);
     HR_wrapper.insert_pair(atom_ij5);
     hamilt::HContainer<double> HR_no_wrapper(&PO);
@@ -628,6 +606,54 @@ TEST_F(HContainerTest, atompair_funcs)
     EXPECT_EQ(HR_no_wrapper.size_R_loop(), 1);
 }
 
+// Test for Wrapper mode in HContainer
+// 1. test constructor of wrapper mode BaseMatrix
+// 2. test constructor of wrapper mode AtomPair
+// 3. test constructor of wrapper mode HContainer
+// 4. test allocate() for wrapper mode
+// 5. test get_nnr() for HContainer
+// 6. test data_access correctnesss for wrapper mode
+TEST_F(HContainerTest, wrapper_mode)
+{
+    EXPECT_EQ(HR->get_nnr(), 36);
+    // test HR_wrapper constructed by HR and vector
+    std::vector<double> hr_data(HR->get_nnr());
+    hamilt::HContainer<double> HR_wrapper(*HR, hr_data.data());
+    EXPECT_EQ(HR_wrapper.get_nnr(), 36);
+    EXPECT_EQ(HR_wrapper.size_atom_pairs(), 9);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(0).get_atom_i(), 0);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(0).get_atom_j(), 0);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(0).get_row_size(), 2);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(0).get_col_size(), 2);
+    EXPECT_EQ(HR_wrapper.get_wrapper(), hr_data.data());
+    EXPECT_EQ(HR_wrapper.get_atom_pair(0).get_pointer(), hr_data.data());
+    for (size_t i = 0; i < hr_data.size(); i++)
+    {
+        hr_data[i] = i;
+    }
+    EXPECT_EQ(HR_wrapper.get_atom_pair(0).get_value(0, 0), 0.0);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(0).get_value(0, 1), 1.0);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(0).get_value(1, 0), 2.0);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(0).get_value(1, 1), 3.0);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(1).get_value(0, 0), 4.0);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(1).get_value(0, 1), 5.0);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(1).get_value(1, 0), 6.0);
+    EXPECT_EQ(HR_wrapper.get_atom_pair(1).get_value(1, 1), 7.0);
+    hamilt::AtomPair<double> atom_ij(HR->get_atom_pair(0), hr_data.data());
+    hamilt::BaseMatrix<double> matrix_test = hamilt::BaseMatrix<double>(atom_ij.get_row_size(), atom_ij.get_col_size(), hr_data.data());
+    EXPECT_EQ(atom_ij.get_value(1, 1), 3.0);
+    EXPECT_EQ(matrix_test.get_value(1, 1), 3.0);
+    HR->allocate(hr_data.data(), false);
+    EXPECT_EQ(HR->get_atom_pair(0).get_value(1, 1), 3.0);
+    EXPECT_EQ(HR->get_atom_pair(1).get_value(1, 1), 7.0);
+    EXPECT_EQ(HR->get_atom_pair(2).get_value(1, 1), 11.0);
+    EXPECT_EQ(HR->get_atom_pair(3).get_value(1, 1), 15.0);
+    HR->allocate(hr_data.data(), true);
+    EXPECT_EQ(HR->get_atom_pair(0).get_value(1, 1), 0.0);
+    EXPECT_EQ(HR->get_atom_pair(1).get_value(1, 1), 0.0);
+    EXPECT_EQ(HR->get_atom_pair(2).get_value(1, 1), 0.0);
+    EXPECT_EQ(HR->get_atom_pair(3).get_value(1, 1), 0.0);
+}
 
 int main(int argc, char** argv)
 {
diff --git a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp
index 3cdfdfa954c..dbc0088ecd9 100644
--- a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp
@@ -59,7 +59,7 @@ TEST_F(HContainerTest, insert_pair)
     // set up a hamilt::AtomPair
     hamilt::AtomPair<std::complex<double>> atom_ij(0, 3);
     atom_ij.set_size(2, 2);
-    atom_ij.allocate(true);
+    atom_ij.allocate(nullptr, true);
     // insert atom_ij into HR
     HR->insert_pair(atom_ij);
     // check if atom_ij is inserted into HR
@@ -83,7 +83,7 @@ TEST_F(HContainerTest, insert_pair)
     // set up another hamilt::AtomPair
     hamilt::AtomPair<std::complex<double>> atom_kl(1, 0);
     atom_kl.set_size(2, 2);
-    atom_kl.allocate(true);
+    atom_kl.allocate(nullptr, true);
     std::complex<double> tmp_array[4] = {1, 2, 3, 4};
     atom_kl.get_HR_values(1, 0, 0).add_array(&tmp_array[0]);
     // insert atom_kl into HR
@@ -285,7 +285,7 @@ TEST_F(HContainerTest, size_atom_pairs)
     // fix to another R
     hamilt::AtomPair<std::complex<double>> atom_ij(0, 1);
     atom_ij.set_size(2, 2);
-    atom_ij.allocate(true);
+    atom_ij.allocate(nullptr, true);
     hamilt::BaseMatrix<std::complex<double>>& tmp = atom_ij.get_HR_values(1, 0, 0);
     std::complex<double> tmp_array[4] = {1, 2, 3, 4};
     tmp.add_array(tmp_array);
@@ -343,7 +343,7 @@ TEST_F(HContainerTest, data)
     // set up a hamilt::AtomPair
     hamilt::AtomPair<std::complex<double>> atom_ij(0, 1);
     atom_ij.set_size(2, 2);
-    atom_ij.allocate(true);
+    atom_ij.allocate(nullptr, true);
     hamilt::BaseMatrix<std::complex<double>>& tmp = atom_ij.get_HR_values(0, 0, 0);
     std::complex<double> tmp_array[4] = {1, 2, 3, 4};
     tmp.add_array(tmp_array);
@@ -374,8 +374,6 @@ TEST_F(HContainerTest, data)
 
 // using TEST_F to test functions in BaseMatrix
 // 1. test constructor with existed data
-// 2. test set_ldc
-// 3. test add_array with memory_type = 2
 // 4. test add_element
 // 5. test get_value
 TEST_F(HContainerTest, basematrix_funcs)
@@ -397,38 +395,6 @@ TEST_F(HContainerTest, basematrix_funcs)
     EXPECT_EQ(BM_copy.get_value(0, 1), std::complex<double>(2));
     EXPECT_EQ(BM_copy.get_value(1, 0), std::complex<double>(3));
     EXPECT_EQ(BM_copy.get_value(1, 1), std::complex<double>(4));
-    // 2. test set_ldc
-    BM.set_ldc(3);
-    EXPECT_EQ(BM.get_value(0, 0), std::complex<double>(1));
-    EXPECT_EQ(BM.get_value(1, 0), std::complex<double>(4));
-    BM.set_ldc(2);
-    EXPECT_EQ(BM.get_value(1, 0), std::complex<double>(3));
-    // 3. test add_array with memory_type = 2
-    std::complex<double> data_ptr2[4] = {5, 6, 7, 8};
-    BM.add_array(&data_ptr2[0]);
-    // check if data is correct
-    EXPECT_EQ(BM.get_value(0, 0), std::complex<double>(6));
-    EXPECT_EQ(data_ptr[0], std::complex<double>(6));
-    EXPECT_EQ(BM.get_value(0, 1), std::complex<double>(8));
-    EXPECT_EQ(BM.get_value(1, 0), std::complex<double>(10));
-    EXPECT_EQ(BM.get_value(1, 1), std::complex<double>(12));
-    // 4. test add_element
-    BM.add_element(0, 0, 1);
-    BM.add_element(0, 1, 2);
-    BM.add_element(1, 0, 3);
-    BM.add_element(1, 1, 4);
-    // check if data is correct
-    EXPECT_EQ(BM.get_value(0, 0), std::complex<double>(7));
-    EXPECT_EQ(BM.get_value(0, 1), std::complex<double>(10));
-    EXPECT_EQ(BM.get_value(1, 0), std::complex<double>(13));
-    EXPECT_EQ(BM.get_value(1, 1), std::complex<double>(16));
-    // 5. test get_value
-    hamilt::BaseMatrix<std::complex<double>> BM2(2, 2);
-    BM2.allocate(true);
-    BM2.add_element(0, 0, 1);
-    BM2.add_element(0, 1, 2);
-    EXPECT_EQ(BM2.get_value(0, 0), std::complex<double>(1));
-    EXPECT_EQ(BM2.get_value(0, 1), std::complex<double>(2));
 } 
 
 // using TEST_F to test functions in AtomPair
@@ -486,7 +452,7 @@ TEST_F(HContainerTest, atompair_funcs)
         for(int atom_j = 0; atom_j<2; ++atom_j)
         {
             hamilt::AtomPair<std::complex<double>> tmp(atom_i, atom_j, 0, 0, 0, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, nullptr);
-            tmp.allocate(false);
+            tmp.allocate(nullptr, false);
             std::complex<double>* tmp_data = tmp.get_HR_values(0, 0, 0).get_pointer();
             for(int i=0;i<4;++i)
             {
@@ -542,7 +508,13 @@ TEST_F(HContainerTest, atompair_funcs)
     HR.unfix_R();
     // check hr_array and hr_array2 are correct
     std::complex<double> correct1;
-    std::complex<double> correct_array[16] = {1, 2, 3, 4, 5, 6, 7, 8, 3, 4, 5, 6, 7, 8, 9, 10};
+    std::complex<double> correct_array[16] = {
+        1, 2, 3, 4, 
+        5, 6, 7, 8, 
+        3, 4, 5, 6, 
+        7, 8, 9, 10};
+    std::complex<double> test_array[16] = {1, 2, 5, 6, 3, 4, 7, 8, 3, 4, 7, 8, 5, 6, 9, 10};
+
     for(int i=0;i<4;++i)
     {
         for(int j=0;j<4;++j)
@@ -553,31 +525,36 @@ TEST_F(HContainerTest, atompair_funcs)
     }
 
     // construct AtomPair from existed matrix
-    hamilt::AtomPair<std::complex<double>> atom_ij4(0, 0, &PO, correct_array);
+    hamilt::AtomPair<std::complex<double>> atom_ij4(0, 0, &PO, test_array);
     EXPECT_EQ(atom_ij4.get_value(0, 0), correct_array[0]);
     EXPECT_EQ(atom_ij4.get_value(1, 1), correct_array[5]);
     EXPECT_EQ(atom_ij4.get_value(0), correct_array[0]);
-    hamilt::AtomPair<std::complex<double>> atom_ij5(0, 1, 1, 1, 1, &PO, correct_array);
-    hamilt::AtomPair<std::complex<double>> atom_ij6(1, 0, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, correct_array);
-    hamilt::AtomPair<std::complex<double>> atom_ij7(1, 1, 1, 1, 1, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, correct_array);
+    hamilt::AtomPair<std::complex<double>> atom_ij5(0, 1, 1, 1, 1, &PO, &test_array[4]);
+    hamilt::AtomPair<std::complex<double>> atom_ij6(1, 0, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, &test_array[8]);
+    hamilt::AtomPair<std::complex<double>> atom_ij7(1, 1, 1, 1, 1, PO.atom_begin_row.data(), PO.atom_begin_col.data(), 2, &test_array[12]);
     // get_matrix_value will use global2local_row and global2local_col in Parallel_Orbitals
     // so we need to set them
     std::ofstream ofs("test_hcontainer_complex.log");
     PO.set_global2local(4, 4, false, ofs);
-    for(int i=0;i<4;++i)
-    {
-        for(int j=0;j<4;++j)
+    auto checkdata = [&](hamilt::AtomPair<std::complex<double>>& ap_in) {
+        auto data_ij4 = ap_in.get_matrix_values();
+        int* tmp_index = std::get<0>(data_ij4).data();
+        std::complex<double>* tmp_data = std::get<1>(data_ij4);
+        double sum_error = 0.0;
+        for(int irow = tmp_index[0]; irow < tmp_index[0] + tmp_index[1]; ++irow)
         {
-            correct1 = correct_array[i*4+j];
-            hamilt::AtomPair<std::complex<double>>* atom_ij_ptr; 
-            if(i/2==0 && j/2==0){ atom_ij_ptr = &atom_ij4; }
-            else if(i/2==0 && j/2==1){ atom_ij_ptr = &atom_ij5; }
-            else if(i/2==1 && j/2==0){ atom_ij_ptr = &atom_ij6; }
-            else if(i/2==1 && j/2==1){ atom_ij_ptr = &atom_ij7; }
-            else{ atom_ij_ptr = nullptr; }
-            EXPECT_EQ(atom_ij_ptr->get_matrix_value(i, j), correct1);
+            for(int icol = tmp_index[2]; icol < tmp_index[2] + tmp_index[3]; ++icol)
+            {
+                sum_error += std::abs((*tmp_data).real() - correct_array[irow*4+icol].real());
+                tmp_data++;
+            }
         }
-    }
+        return sum_error;
+    };
+    EXPECT_EQ(checkdata(atom_ij4), 0.0);
+    EXPECT_EQ(checkdata(atom_ij5), 0.0);
+    EXPECT_EQ(checkdata(atom_ij6), 0.0);
+    EXPECT_EQ(checkdata(atom_ij7), 0.0);
 }
 
 
diff --git a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_output.cpp b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_output.cpp
index de3a27c7c67..5e37a17bcd2 100644
--- a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_output.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_output.cpp
@@ -91,11 +91,11 @@ TEST_F(OutputHContainerTest, Write)
     // 5 0 7 0
     // 3 0 5 6
     // 7 8 0 10
-    hamilt::AtomPair<double> ap1(0, 1, 0, 1, 1, &ParaV, correct_array);
-    hamilt::AtomPair<double> ap2(1, 1, 0, 0, 0, &ParaV, correct_array1);
+    double test_data[8] = {0, 4, 7, 0, 5, 6, 0, 10};
+    hamilt::AtomPair<double> ap1(0, 1, 0, 1, 1, &ParaV, &test_data[0]);
+    hamilt::AtomPair<double> ap2(1, 1, 0, 0, 0, &ParaV, &test_data[4]);
     HR.insert_pair(ap1);
     HR.insert_pair(ap2);
-    HR.allocate(true);
     for (int ir = 0; ir < HR.size_R_loop(); ++ir)
     {
         int rx, ry, rz;
@@ -111,10 +111,10 @@ TEST_F(OutputHContainerTest, Write)
                 EXPECT_DOUBLE_EQ(tmp_ap.get_value(0, 1), 4);
                 EXPECT_DOUBLE_EQ(tmp_ap.get_value(1, 0), 7);
                 EXPECT_DOUBLE_EQ(tmp_ap.get_value(1, 1), 0);
-                EXPECT_DOUBLE_EQ(tmp_ap.get_matrix_value(0, 2), 0);
-                EXPECT_DOUBLE_EQ(tmp_ap.get_matrix_value(0, 3), 4);
-                EXPECT_DOUBLE_EQ(tmp_ap.get_matrix_value(1, 2), 7);
-                EXPECT_DOUBLE_EQ(tmp_ap.get_matrix_value(1, 3), 0);
+                EXPECT_DOUBLE_EQ(std::get<0>(tmp_ap.get_matrix_values())[0], 0);
+                EXPECT_DOUBLE_EQ(std::get<0>(tmp_ap.get_matrix_values())[1], 2);
+                EXPECT_DOUBLE_EQ(std::get<0>(tmp_ap.get_matrix_values())[2], 2);
+                EXPECT_DOUBLE_EQ(std::get<0>(tmp_ap.get_matrix_values())[3], 2);
             }
             else if (rx == 0 && ry == 0 && rz == 0)
             {
@@ -122,10 +122,10 @@ TEST_F(OutputHContainerTest, Write)
                 EXPECT_DOUBLE_EQ(tmp_ap.get_value(0, 1), 6);
                 EXPECT_DOUBLE_EQ(tmp_ap.get_value(1, 0), 0);
                 EXPECT_DOUBLE_EQ(tmp_ap.get_value(1, 1), 10);
-                EXPECT_DOUBLE_EQ(tmp_ap.get_matrix_value(2, 2), 5);
-                EXPECT_DOUBLE_EQ(tmp_ap.get_matrix_value(2, 3), 6);
-                EXPECT_DOUBLE_EQ(tmp_ap.get_matrix_value(3, 2), 0);
-                EXPECT_DOUBLE_EQ(tmp_ap.get_matrix_value(3, 3), 10);
+                EXPECT_DOUBLE_EQ(std::get<0>(tmp_ap.get_matrix_values())[0], 2);
+                EXPECT_DOUBLE_EQ(std::get<0>(tmp_ap.get_matrix_values())[1], 2);
+                EXPECT_DOUBLE_EQ(std::get<0>(tmp_ap.get_matrix_values())[2], 2);
+                EXPECT_DOUBLE_EQ(std::get<0>(tmp_ap.get_matrix_values())[3], 2);
             }
         }
         HR.unfix_R();
diff --git a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_readCSR.cpp b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_readCSR.cpp
index 76660803a7d..8a4affd1f60 100644
--- a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_readCSR.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_readCSR.cpp
@@ -98,7 +98,7 @@ TEST_F(ReadHContainerTest, ReadAndOutputHContainer)
                 int end_col = paraV.atom_begin_col[jat + 1];
                 int numberofCol = end_col - begin_col;
                 hamilt::BaseMatrix<double> tmp_matrix(numberofRow, numberofCol);
-                tmp_matrix.allocate(true);
+                tmp_matrix.allocate(nullptr, true);
                 int nnz = 0;
                 for (const auto& element: sparse_matrix.getElements())
                 {
@@ -114,7 +114,7 @@ TEST_F(ReadHContainerTest, ReadAndOutputHContainer)
                 if (nnz != 0)
                 {
                     auto tmp_ap = hamilt::AtomPair<double>(iat, jat, RCoord[0], RCoord[1], RCoord[2], &paraV);
-                    tmp_ap.allocate(true);
+                    tmp_ap.allocate(nullptr, true);
                     tmp_ap.convert_add(tmp_matrix, RCoord[0], RCoord[1], RCoord[2]);
                     SR.insert_pair(tmp_ap);
                 }
diff --git a/source/module_hamilt_lcao/module_hcontainer/test/test_transfer.cpp b/source/module_hamilt_lcao/module_hcontainer/test/test_transfer.cpp
index 714ed6467ff..8e07d2942e4 100644
--- a/source/module_hamilt_lcao/module_hcontainer/test/test_transfer.cpp
+++ b/source/module_hamilt_lcao/module_hcontainer/test/test_transfer.cpp
@@ -233,6 +233,42 @@ TEST_F(TransferTest, paraToSerial)
     end_time = std::chrono::high_resolution_clock::now();
     std::chrono::duration<double> elapsed_time2 = std::chrono::duration_cast<std::chrono::duration<double>>(end_time - start_time);
     std::cout <<"rank: "<<my_rank<< " HR init time: " << elapsed_time0.count()<<" transfer_p2s time: "<<elapsed_time1.count()<<" check time: "<<elapsed_time2.count()<<" seconds." << std::endl;
+
+    // now assume the HR_serial is empty and do the same test
+    Parallel_Orbitals serialV;
+    std::ofstream ofs("test.log");
+    serialV.set_global2local(test_size*test_nw, test_size*test_nw, false, ofs);
+    serialV.set_atomic_trace(ucell.get_iat2iwt(), test_size, test_size*test_nw);
+    hamilt::HContainer<double>* HR_serial2;
+    if(my_rank == 0)
+    {
+        HR_serial2 = new hamilt::HContainer<double>(&serialV);
+    }
+    hamilt::gatherParallels(*HR_para, HR_serial2, 0);
+    if(my_rank == 0)
+    {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+        for(int i = 0; i < HR_serial2->size_atom_pairs(); i++)
+        {
+            hamilt::AtomPair<double>& atom_pair = HR_serial2->get_atom_pair(i);
+            int atom_i = atom_pair.get_atom_i();
+            int atom_j = atom_pair.get_atom_j();
+            //lambda function to calculate value of array: (atom_i*test_size+atom_j+k)*test_nw + l
+            auto value = [&](int k, int l) -> double {return (((atom_i*test_nw+k)*test_size+atom_j)*test_nw + l);};
+            double* data = atom_pair.get_pointer(0);
+            for(int k = 0; k < test_nw; k++)
+            {
+                for(int l = 0; l < test_nw; l++)
+                {
+                    EXPECT_NEAR(*data , value(k, l), 1e-10);
+                    ++data;
+                }
+            }
+        }
+        delete HR_serial2;
+    }
 #endif
 }