diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index 2e190877d92..768e6238a4a 100644
--- a/source/Makefile.Objects
+++ b/source/Makefile.Objects
@@ -393,6 +393,7 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
       write_dm_sparse.o\
 
 OBJS_LCAO=DM_gamma.o\
+      DM_gamma_2d_to_grid.o\
       DM_k.o\
       evolve_elec.o\
       evolve_psi.o\
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt
index 5d5a9ad7da1..a05ffa79471 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt
@@ -27,6 +27,7 @@ if(ENABLE_LCAO)
         center2_orb-orb21.cpp
         center2_orb-orb22.cpp
         DM_gamma.cpp
+        DM_gamma_2d_to_grid.cpp
         DM_k.cpp
         local_orbital_charge.cpp
         local_orbital_wfc.cpp
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma.cpp
index 715c2427e6c..56a1b194306 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma.cpp
@@ -1,208 +1,11 @@
 #include "local_orbital_charge.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
-#include "module_base/blas_connector.h"
 #include "module_io/read_wfc_nao.h"
 #include "module_base/parallel_reduce.h"
 #include "module_base/parallel_common.h"
 #include "module_base/memory.h"
 #include "module_base/timer.h"
 
-extern "C"
-{
-    void Cblacs_gridinfo(int icontxt, int* nprow, int *npcol, int *myprow, int *mypcol);
-    void Cblacs_pinfo(int *myid, int *nprocs);
-    void Cblacs_pcoord(int icontxt, int pnum, int *prow, int *pcol);
-    int Cblacs_pnum(int icontxt, int prow, int pcol);
-}
-
-#ifdef __MPI
-int Local_Orbital_Charge::setAlltoallvParameter(MPI_Comm comm_2D, int blacs_ctxt, int nblk)
-{
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"enter setAlltoallvParameter, nblk", nblk);
-    ModuleBase::timer::tick("LOC","Alltoall");
-    // setup blacs parameters
-    int nprows=0;	
-	int npcols=0;
-	int nprocs=0;
-    int myprow=0;
-	int mypcol=0;
-	int myproc=0;
-
-    Cblacs_gridinfo(blacs_ctxt, &nprows, &npcols, &myprow, &mypcol);
-
-    Cblacs_pinfo(&myproc, &nprocs);
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nprocs",nprocs);
-
-    size_t memory_sender = 0;
-    size_t memory_receiver = 0;
-    size_t memory_other = 0;
-    // init data arrays
-    delete[] sender_size_process;
-    sender_size_process=new int[nprocs];
-    delete[] sender_displacement_process;
-    sender_displacement_process=new int[nprocs];
-    //GlobalV::ofs_running << "checkpoint 2" << std::endl;
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"lgd_now",lgd_now);
-    memory_sender += nprocs * 2 * sizeof(int);
-    
-    receiver_size=lgd_now*lgd_now;
-    delete[] receiver_size_process;
-    receiver_size_process=new int[nprocs];
-    delete[] receiver_displacement_process;
-    receiver_displacement_process=new int[nprocs];
-    delete[] receiver_local_index;
-    receiver_local_index=new int[receiver_size];
-    delete[] receiver_buffer;
-    receiver_buffer=new double[receiver_size];
-
-    memory_receiver += nprocs * 2 * sizeof(int) + receiver_size * (sizeof(int) + sizeof(double));
-    
-    int *trace_2D_row=new int[lgd_now];
-    int *trace_2D_col=new int[lgd_now];
-    int *trace_2D_prow=new int[lgd_now];
-    int *trace_2D_pcol=new int[lgd_now];
-    //int *trace_global=new int[lgd_now];
-
-    int *nRow_in_proc=new int[nprows];
-    int *nCol_in_proc=new int[npcols];
-
-    memory_other += (lgd_now * 4 + nprows + npcols) * sizeof(int);
-
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nprows",nprows);
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"npcols",npcols);
-
-    for(int i=0; i<nprows; ++i)
-    {
-        nRow_in_proc[i]=0;
-    }
-    for(int i=0; i<npcols; ++i)
-    {
-        nCol_in_proc[i]=0;
-    }
-
-    // count the number of elements to be received from each process
-    for(int iGlobal=0; iGlobal<GlobalV::NLOCAL; ++iGlobal)
-    {
-        int iLocalGrid=this->LOWF->gridt->trace_lo[iGlobal];
-        if(iLocalGrid>=0)
-        {
-            //trace_global[iLocalGrid]=iGlobal;
-            int p;
-            trace_2D_row[iLocalGrid]=Local_Orbital_wfc::localIndex(iGlobal, nblk, nprows, p);
-            trace_2D_prow[iLocalGrid]=p;
-            nRow_in_proc[trace_2D_prow[iLocalGrid]]++;
-            trace_2D_col[iLocalGrid]=Local_Orbital_wfc::localIndex(iGlobal, nblk, npcols, p);
-            trace_2D_pcol[iLocalGrid]=p;
-            nCol_in_proc[trace_2D_pcol[iLocalGrid]]++;
-        }
-    }
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"NLOCAL",GlobalV::NLOCAL);
-    receiver_displacement_process[0]=0;
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"receiver_displacement_process[0]",receiver_displacement_process[0]);
-    for(int pnum=0; pnum<nprocs; ++pnum)
-    {
-        int prow, pcol;
-        Cblacs_pcoord(blacs_ctxt, pnum, &prow, &pcol);
-        receiver_size_process[pnum]=nRow_in_proc[prow]*nCol_in_proc[pcol];
-
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"pnum",pnum);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"prow",prow);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"pcol",pcol);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nRow_in_proc",nRow_in_proc[prow]);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nCol_in_proc",nCol_in_proc[pcol]);
-
-        if(pnum>0)
-        {
-            receiver_displacement_process[pnum]=receiver_displacement_process[pnum-1]+receiver_size_process[pnum-1];
-        }
-    }
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last receiver_size_process",receiver_size_process[nprocs-1]);
-    
-    // build the index to be received
-    int* pos=new int[nprocs];
-    int *receiver_2D_index=new int[receiver_size];
-    for(int i=0; i<nprocs; ++i)
-    {
-        pos[i]=receiver_displacement_process[i];
-    }
-    memory_other += nprocs * sizeof(int);
-    memory_receiver += receiver_size * sizeof(int);
-
-    for(int i=0; i<lgd_now; ++i)
-    {
-        int src_row=trace_2D_row[i];
-        int src_prow=trace_2D_prow[i];
-        for(int j=0; j<lgd_now; ++j)
-        {
-            int src_col=trace_2D_col[j];
-            int src_idx=src_row*GlobalV::NLOCAL+src_col; // leanding dimension is set to GlobalV::NLOCAL for all processes
-
-            int src_pcol=trace_2D_pcol[j];
-            int src_proc=Cblacs_pnum(blacs_ctxt, src_prow, src_pcol);
-
-            receiver_2D_index[pos[src_proc]]=src_idx;
-            receiver_local_index[pos[src_proc]]=i*lgd_now+j;
-            ++pos[src_proc];
-        }
-    }
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last receiver_2D_index",receiver_2D_index[lgd_now*lgd_now-1]);
-    delete[] pos;
-    delete[] trace_2D_row;
-    delete[] trace_2D_col;
-    delete[] trace_2D_prow;
-    delete[] trace_2D_pcol;
-    //delete[] trace_global;
-    delete[] nRow_in_proc;
-    delete[] nCol_in_proc;
-    
-    // send number of elements to be sent via MPI_Alltoall
-    MPI_Alltoall(receiver_size_process, 1, MPI_INT,
-                 sender_size_process, 1, MPI_INT, comm_2D);
-    
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last sender_size_process",sender_size_process[nprocs-1]);
-    // setup sender buffer
-    sender_size=sender_size_process[0];
-    sender_displacement_process[0]=0;
-    for(int i=1; i<nprocs; ++i)
-    {
-        sender_size+=sender_size_process[i];
-        sender_displacement_process[i]=sender_displacement_process[i-1]+sender_size_process[i-1];
-    }
-    
-    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"sender_size",sender_size);
-    delete[] sender_2D_index;
-    sender_2D_index=new int[sender_size];
-    delete[] sender_buffer;
-    sender_buffer=new double[sender_size];
-    memory_sender += sender_size * (sizeof(int) * sizeof(double));
-
-    // send the index of the elements to be received via MPI_Alltoall
-    MPI_Alltoallv(receiver_2D_index, receiver_size_process, receiver_displacement_process, MPI_INT,
-                  sender_2D_index, sender_size_process, sender_displacement_process, MPI_INT, comm_2D);
-
-
-    GlobalV::ofs_running << "receiver_size is " << receiver_size << " ; receiver_size of each process is:\n";
-    for(int i=0; i<nprocs; ++i)
-    {
-        GlobalV::ofs_running<<receiver_size_process[i]<<" ";
-    }
-    GlobalV::ofs_running<<std::endl;
-    GlobalV::ofs_running<<"sender_size is "<<sender_size<<" ; sender_size of each process is:\n";
-    for(int i=0; i<nprocs; ++i)
-    {
-        GlobalV::ofs_running<<sender_size_process[i]<<" ";
-    }
-    GlobalV::ofs_running << std::endl;
-        
-        // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last sender_2D_index",sender_2D_index[lgd_now*lgd_now-1]);
-    delete[] receiver_2D_index;
-    ModuleBase::Memory::record("LOC::A2A_receiv", memory_receiver);
-    ModuleBase::Memory::record("LOC::A2A_sender", memory_sender);
-    ModuleBase::Memory::record("LOC::A2A_other", memory_other);
-    ModuleBase::timer::tick("LOC","Alltoall");
-    return 0;
-}
-#endif
 
 // allocate density kernel may change once the ion
 // positions change
@@ -268,7 +71,7 @@ void Local_Orbital_Charge::allocate_gamma(
     }
     
 #ifdef __MPI
-    setAlltoallvParameter(this->ParaV->comm_2D, this->ParaV->blacs_ctxt, this->ParaV->nb);
+    this->dm2g.setAlltoallvParameter(this->ParaV->comm_2D, GlobalV::NLOCAL, this->ParaV->blacs_ctxt, this->ParaV->nb, this->lgd_now, this->LOWF->gridt->trace_lo);
 #endif
 
 	// Peize Lin test 2019-01-16
@@ -343,124 +146,6 @@ void Local_Orbital_Charge::cal_dk_gamma_from_2D_pub(void)
 {
     ModuleBase::TITLE("Local_Orbital_Charge","cal_dk_gamma_from_2D_pub");
 
-	cal_dk_gamma_from_2D();
+    this->dm2g.cal_dk_gamma_from_2D(this->dm_gamma, this->DM, GlobalV::NSPIN, GlobalV::NLOCAL, this->lgd_now, GlobalV::ofs_running);
 }
-// calculate the grid distributed DM matrix from 2D block-cyclic distributed DM matrix
-// transform dm_gamma[is].c to this->DM[is]
-void Local_Orbital_Charge::cal_dk_gamma_from_2D(void)
-{
-    ModuleBase::timer::tick("LCAO_Charge","dm_2dTOgrid");
-#ifdef __DEBUG
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"cal_dk_gamma_from_2D, NSPIN", GlobalV::NSPIN);
-#endif
-
-    for(int is=0; is<GlobalV::NSPIN; ++is)
-    {
-        // int myid;
-        // MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-        // if(myid==0)
-        // {
-        //     GlobalV::ofs_running<<"DM[0][0:1][0:1] before send:"<<std::endl;
-        //     GlobalV::ofs_running<<"DM(0,0)"<<wfc_dm_2d.dm_gamma[is](0,0)<<" ";
-        //     GlobalV::ofs_running<<"DM(0,1)"<<wfc_dm_2d.dm_gamma[is](1,0)<<std::endl;
-        //     GlobalV::ofs_running<<"DM(1,0)"<<wfc_dm_2d.dm_gamma[is](0,1)<<" ";
-        //     GlobalV::ofs_running<<"DM(1,1)"<<wfc_dm_2d.dm_gamma[is](1,1)<<std::endl;
-        // }
-        /*GlobalV::ofs_running<<"2D block parameters:\n"<<"nblk: "<<this->ParaV->nb<<std::endl;
-        GlobalV::ofs_running<<"DM in 2D format:\n_________________________________________\n";
-        for(int i=0; i<this->dm_gamma[is].nr; ++i)
-        {
-            for(int j=0; j<this->dm_gamma[is].nc; ++j)
-            {
-                GlobalV::ofs_running<<this->dm_gamma[is](i,j)<<" ";
-            }
-            GlobalV::ofs_running<<std::endl;
-        }
-        GlobalV::ofs_running<<"=========================================\n";*/
-
-        // put data from dm_gamma[is] to sender index
-        int nNONZERO=0;
-        for(int i=0; i<sender_size; ++i)
-        {
-            const int idx=sender_2D_index[i];
-            const int icol=idx%GlobalV::NLOCAL;
-            const int irow=(idx-icol)/GlobalV::NLOCAL;
-            // sender_buffer[i]=wfc_dm_2d.dm_gamma[is](irow,icol);
-            sender_buffer[i]=this->dm_gamma[is](icol,irow); // sender_buffer is clomun major, 
-                                                                // so the row and column index should be switched
-            if(sender_buffer[i]!=0) ++nNONZERO;
-        }
 
-#ifdef __DEBUG
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"number of non-zero elements in sender_buffer",nNONZERO);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"sender_size",sender_size);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last sender_buffer",sender_buffer[sender_size-1]);
-#endif
-
-        // transform data via MPI_Alltoallv
-        #ifdef __MPI
-        MPI_Alltoallv(sender_buffer, sender_size_process, sender_displacement_process, MPI_DOUBLE,
-                      receiver_buffer, receiver_size_process, receiver_displacement_process, MPI_DOUBLE, this->ParaV->comm_2D);
-        #endif
-        // put data from receiver buffer to this->DM[is]
-        nNONZERO=0;
-        // init DM[is]
-        /*for(int i=0; i<lgd_now; ++i)
-        {
-            for(int j=0; j<lgd_now; ++j)
-            {
-                DM[is][i][j]=0;
-            }
-        }*/
-        for(int i=0; i<receiver_size; ++i)
-        {
-            const int idx=receiver_local_index[i];
-            const int icol=idx%lgd_now;
-            const int irow=(idx-icol)/lgd_now;
-            DM[is][irow][icol]=receiver_buffer[i];
-            //DM[is][icol][irow]=receiver_buffer[i];
-            if(receiver_buffer[i]!=0) ++nNONZERO;
-        }
-
-#ifdef __DEBUG
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"number of non-zero elements in receiver_buffer",nNONZERO);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"receiver_size",receiver_size);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last receiver_buffer",receiver_buffer[receiver_size-1]);
-#endif
-        // GlobalV::ofs_running<<"DM[0][0:1][0:1] after receiver:"<<std::endl;
-        // int idx0=GlobalC::GridT.trace_lo[0];
-        // int idx1=GlobalC::GridT.trace_lo[1];
-        // if(idx0>=0)
-        // {
-        //     GlobalV::ofs_running<<"DM(0,0)"<<DM[0][idx0][idx0]<<" ";
-        // }
-        // if(idx0>=0 && idx1>=0)
-        // {
-        //     GlobalV::ofs_running<<"DM(0,1)"<<DM[0][idx0][idx1]<<std::endl;
-        //     GlobalV::ofs_running<<"DM(1,0)"<<DM[0][idx1][idx0]<<" ";
-        // }
-        // if(idx1>=0)
-        // {
-        //     GlobalV::ofs_running<<"DM(1,1)"<<DM[0][idx1][idx1]<<std::endl;
-        // }
-        //GlobalV::ofs_running<<DM[0][0][0]<<" "<<DM[0][0][1]<<std::endl;
-        //GlobalV::ofs_running<<DM[0][1][0]<<" "<<DM[0][1][1]<<std::endl;
-        /*GlobalV::ofs_running<<"DM in local grid:\n_________________________________________\n";
-        for(int i=0; i<GlobalV::NLOCAL; ++i)
-        {
-            int ii=GlobalC::GridT.trace_lo[i];
-            if(ii < 0) continue;
-            for(int j=0; j<GlobalV::NLOCAL; ++j)
-            {
-                int jj=GlobalC::GridT.trace_lo[j];
-                if(jj<0) continue;
-                GlobalV::ofs_running<<DM[is][ii][jj]<<" ";
-            }
-            GlobalV::ofs_running<<std::endl;
-        }
-        GlobalV::ofs_running<<"=========================================\n";*/
-
-    }
-    ModuleBase::timer::tick("LCAO_Charge","dm_2dTOgrid");
-	return;
-}
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.cpp
new file mode 100644
index 00000000000..24977009370
--- /dev/null
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.cpp
@@ -0,0 +1,290 @@
+#include "DM_gamma_2d_to_grid.h"
+#include "module_base/global_function.h"
+#include "module_base/memory.h"
+#include "module_base/timer.h"
+#ifdef __MPI
+#include "module_base/blacs_connector.h"
+#endif
+
+DMgamma_2dtoGrid::DMgamma_2dtoGrid()
+{
+    sender_2D_index = nullptr;
+    sender_size_process = nullptr;
+    sender_displacement_process = nullptr;
+    sender_buffer = nullptr;
+
+    receiver_local_index = nullptr;
+    receiver_size_process = nullptr;
+    receiver_displacement_process = nullptr;
+    receiver_buffer = nullptr;
+}
+DMgamma_2dtoGrid::~DMgamma_2dtoGrid()
+{
+    delete[] sender_2D_index;
+    delete[] sender_size_process;
+    delete[] sender_displacement_process;
+    delete[] sender_buffer;
+
+    delete[] receiver_local_index;
+    delete[] receiver_size_process;
+    delete[] receiver_displacement_process;
+    delete[] receiver_buffer;
+}
+#ifdef __MPI
+int DMgamma_2dtoGrid::setAlltoallvParameter(MPI_Comm comm_2D, int nbasis, int blacs_ctxt, int nblk, const int& loc_grid_dim, const int* global2local_grid)
+{
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "enter setAlltoallvParameter, nblk", nblk);
+    ModuleBase::timer::tick("LOC", "Alltoall");
+    this->comm_2D = comm_2D;
+    // setup blacs parameters
+    int nprows = 0;
+    int npcols = 0;
+    int nprocs = 0;
+    int myprow = 0;
+    int mypcol = 0;
+    int myproc = 0;
+
+    Cblacs_gridinfo(blacs_ctxt, &nprows, &npcols, &myprow, &mypcol);
+
+    Cblacs_pinfo(&myproc, &nprocs);
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nprocs",nprocs);
+
+    size_t memory_sender = 0;
+    size_t memory_receiver = 0;
+    size_t memory_other = 0;
+    // init data arrays
+    delete[] sender_size_process;
+    sender_size_process = new int[nprocs];
+    delete[] sender_displacement_process;
+    sender_displacement_process = new int[nprocs];
+    //GlobalV::ofs_running << "checkpoint 2" << std::endl;
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"lgd_now",lgd_now);
+    memory_sender += nprocs * 2 * sizeof(int);
+
+    receiver_size = loc_grid_dim * loc_grid_dim;
+    delete[] receiver_size_process;
+    receiver_size_process = new int[nprocs];
+    delete[] receiver_displacement_process;
+    receiver_displacement_process = new int[nprocs];
+    delete[] receiver_local_index;
+    receiver_local_index = new int[receiver_size];
+    delete[] receiver_buffer;
+    receiver_buffer = new double[receiver_size];
+
+    memory_receiver += nprocs * 2 * sizeof(int) + receiver_size * (sizeof(int) + sizeof(double));
+
+    int* trace_2D_row = new int[loc_grid_dim];
+    int* trace_2D_col = new int[loc_grid_dim];
+    int* trace_2D_prow = new int[loc_grid_dim];
+    int* trace_2D_pcol = new int[loc_grid_dim];
+
+    int* nRow_in_proc = new int[nprows];
+    int* nCol_in_proc = new int[npcols];
+
+    memory_other += (loc_grid_dim * 4 + nprows + npcols) * sizeof(int);
+
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nprows",nprows);
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"npcols",npcols);
+
+    for (int i = 0; i < nprows; ++i)
+    {
+        nRow_in_proc[i] = 0;
+    }
+    for (int i = 0; i < npcols; ++i)
+    {
+        nCol_in_proc[i] = 0;
+    }
+
+    // count the number of elements to be received from each process
+    auto localIndex = [](int globalindex, int nblk, int nprocs, int& myproc) -> int
+        {
+            myproc = int((globalindex % (nblk * nprocs)) / nblk);
+            return int(globalindex / (nblk * nprocs)) * nblk + globalindex % nblk;
+        };
+    for (int iGlobal = 0; iGlobal < nbasis; ++iGlobal)
+    {
+        int iLocalGrid = global2local_grid[iGlobal];
+        if (iLocalGrid >= 0)
+        {
+            //trace_global[iLocalGrid]=iGlobal;
+            int p;
+            trace_2D_row[iLocalGrid] = localIndex(iGlobal, nblk, nprows, p);
+            trace_2D_prow[iLocalGrid] = p;
+            nRow_in_proc[trace_2D_prow[iLocalGrid]]++;
+            trace_2D_col[iLocalGrid] = localIndex(iGlobal, nblk, npcols, p);
+            trace_2D_pcol[iLocalGrid] = p;
+            nCol_in_proc[trace_2D_pcol[iLocalGrid]]++;
+        }
+    }
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"NLOCAL",nbasis);
+    receiver_displacement_process[0] = 0;
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"receiver_displacement_process[0]",receiver_displacement_process[0]);
+    for (int pnum = 0; pnum < nprocs; ++pnum)
+    {
+        int prow, pcol;
+        Cblacs_pcoord(blacs_ctxt, pnum, &prow, &pcol);
+        receiver_size_process[pnum] = nRow_in_proc[prow] * nCol_in_proc[pcol];
+
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "pnum", pnum);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "prow", prow);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "pcol", pcol);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nRow_in_proc", nRow_in_proc[prow]);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nCol_in_proc", nCol_in_proc[pcol]);
+
+        if (pnum > 0)
+        {
+            receiver_displacement_process[pnum] = receiver_displacement_process[pnum - 1] + receiver_size_process[pnum - 1];
+        }
+    }
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last receiver_size_process",receiver_size_process[nprocs-1]);
+
+    // build the index to be received
+    int* pos = new int[nprocs];
+    int* receiver_2D_index = new int[receiver_size];
+    for (int i = 0; i < nprocs; ++i)
+    {
+        pos[i] = receiver_displacement_process[i];
+    }
+    memory_other += nprocs * sizeof(int);
+    memory_receiver += receiver_size * sizeof(int);
+
+    for (int i = 0; i < loc_grid_dim; ++i)
+    {
+        int src_row = trace_2D_row[i];
+        int src_prow = trace_2D_prow[i];
+        for (int j = 0; j < loc_grid_dim; ++j)
+        {
+            int src_col = trace_2D_col[j];
+            int src_idx = src_row * nbasis + src_col; // leanding dimension is set to nbasis for all processes
+
+            int src_pcol = trace_2D_pcol[j];
+            int src_proc = Cblacs_pnum(blacs_ctxt, src_prow, src_pcol);
+
+            receiver_2D_index[pos[src_proc]] = src_idx;
+            receiver_local_index[pos[src_proc]] = i * loc_grid_dim + j;
+            ++pos[src_proc];
+        }
+    }
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last receiver_2D_index",receiver_2D_index[lgd_now*lgd_now-1]);
+    delete[] pos;
+    delete[] trace_2D_row;
+    delete[] trace_2D_col;
+    delete[] trace_2D_prow;
+    delete[] trace_2D_pcol;
+    //delete[] trace_global;
+    delete[] nRow_in_proc;
+    delete[] nCol_in_proc;
+
+    // send number of elements to be sent via MPI_Alltoall
+    MPI_Alltoall(receiver_size_process, 1, MPI_INT,
+        sender_size_process, 1, MPI_INT, comm_2D);
+
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last sender_size_process",sender_size_process[nprocs-1]);
+    // setup sender buffer
+    sender_size = sender_size_process[0];
+    sender_displacement_process[0] = 0;
+    for (int i = 1; i < nprocs; ++i)
+    {
+        sender_size += sender_size_process[i];
+        sender_displacement_process[i] = sender_displacement_process[i - 1] + sender_size_process[i - 1];
+    }
+
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"sender_size",sender_size);
+    delete[] sender_2D_index;
+    sender_2D_index = new int[sender_size];
+    delete[] sender_buffer;
+    sender_buffer = new double[sender_size];
+    memory_sender += sender_size * (sizeof(int) * sizeof(double));
+
+    // send the index of the elements to be received via MPI_Alltoall
+    MPI_Alltoallv(receiver_2D_index, receiver_size_process, receiver_displacement_process, MPI_INT,
+        sender_2D_index, sender_size_process, sender_displacement_process, MPI_INT, comm_2D);
+
+
+    GlobalV::ofs_running << "receiver_size is " << receiver_size << " ; receiver_size of each process is:\n";
+    for (int i = 0; i < nprocs; ++i)
+    {
+        GlobalV::ofs_running << receiver_size_process[i] << " ";
+    }
+    GlobalV::ofs_running << std::endl;
+    GlobalV::ofs_running << "sender_size is " << sender_size << " ; sender_size of each process is:\n";
+    for (int i = 0; i < nprocs; ++i)
+    {
+        GlobalV::ofs_running << sender_size_process[i] << " ";
+    }
+    GlobalV::ofs_running << std::endl;
+
+    // ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"last sender_2D_index",sender_2D_index[lgd_now*lgd_now-1]);
+    delete[] receiver_2D_index;
+    ModuleBase::Memory::record("LOC::A2A_receiv", memory_receiver);
+    ModuleBase::Memory::record("LOC::A2A_sender", memory_sender);
+    ModuleBase::Memory::record("LOC::A2A_other", memory_other);
+    ModuleBase::timer::tick("LOC", "Alltoall");
+    return 0;
+}
+#endif
+
+// calculate the grid distributed DM matrix from 2D block-cyclic distributed DM matrix
+// transform dm_gamma[is].c to this->DM[is]
+void DMgamma_2dtoGrid::cal_dk_gamma_from_2D(
+    std::vector<ModuleBase::matrix>& dm_gamma_2d,
+    double*** dm_gamma_grid,
+    int& nspin,
+    int& nbasis,
+    int& loc_grid_dim,
+    std::ofstream& ofs_running)
+{
+    ModuleBase::timer::tick("LCAO_Charge", "dm_2dTOgrid");
+#ifdef __DEBUG
+    ModuleBase::GlobalFunc::OUT(ofs_running, "cal_dk_gamma_from_2D, NSPIN", NSPIN);
+#endif
+
+    for (int is = 0; is < nspin; ++is)
+    {
+        // put data from dm_gamma[is] to sender index
+        int nNONZERO = 0;
+        for (int i = 0; i < sender_size; ++i)
+        {
+            const int idx = sender_2D_index[i];
+            const int icol = idx % nbasis;
+            const int irow = (idx - icol) / nbasis;
+            // sender_buffer[i]=wfc_dm_2d.dm_gamma[is](irow,icol);
+            sender_buffer[i] = dm_gamma_2d[is](icol, irow); // sender_buffer is clomun major, 
+            // so the row and column index should be switched
+            if (sender_buffer[i] != 0) ++nNONZERO;
+        }
+
+#ifdef __DEBUG
+        ModuleBase::GlobalFunc::OUT(ofs_running, "number of non-zero elements in sender_buffer", nNONZERO);
+        ModuleBase::GlobalFunc::OUT(ofs_running, "sender_size", sender_size);
+        ModuleBase::GlobalFunc::OUT(ofs_running, "last sender_buffer", sender_buffer[sender_size - 1]);
+#endif
+
+        // transform data via MPI_Alltoallv
+#ifdef __MPI
+        MPI_Alltoallv(sender_buffer, sender_size_process, sender_displacement_process, MPI_DOUBLE,
+            receiver_buffer, receiver_size_process, receiver_displacement_process, MPI_DOUBLE, this->comm_2D);
+#endif
+        // put data from receiver buffer to this->DM[is]
+        nNONZERO = 0;
+
+        for (int i = 0; i < receiver_size; ++i)
+        {
+            const int idx = receiver_local_index[i];
+            const int icol = idx % loc_grid_dim;
+            const int irow = (idx - icol) / loc_grid_dim;
+            dm_gamma_grid[is][irow][icol] = receiver_buffer[i];
+            //DM[is][icol][irow]=receiver_buffer[i];
+            if (receiver_buffer[i] != 0) ++nNONZERO;
+        }
+
+#ifdef __DEBUG
+        ModuleBase::GlobalFunc::OUT(ofs_running, "number of non-zero elements in receiver_buffer", nNONZERO);
+        ModuleBase::GlobalFunc::OUT(ofs_running, "receiver_size", receiver_size);
+        ModuleBase::GlobalFunc::OUT(ofs_running, "last receiver_buffer", receiver_buffer[receiver_size - 1]);
+#endif
+
+    }
+    ModuleBase::timer::tick("LCAO_Charge", "dm_2dTOgrid");
+    return;
+}
\ No newline at end of file
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.h b/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.h
new file mode 100644
index 00000000000..a38a1f2de06
--- /dev/null
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.h
@@ -0,0 +1,41 @@
+#include <vector>
+#include "module_base/matrix.h"
+#ifdef __MPI
+#include <mpi.h>
+#endif
+
+
+
+class DMgamma_2dtoGrid
+{
+public:
+    DMgamma_2dtoGrid();
+    ~DMgamma_2dtoGrid();
+#ifdef __MPI
+    int setAlltoallvParameter(MPI_Comm comm_2D, int nbasis, int blacs_ctxt, int nblk, const int& loc_grid_dim, const int* global2local_grid);
+#endif
+    void cal_dk_gamma_from_2D(
+        std::vector<ModuleBase::matrix>& dm_gamma_2d,
+        double*** dm_gamma_grid,
+        int& nspin,
+        int& nbasis,
+        int& loc_grid_dim,
+        std::ofstream& ofs_running);
+private:
+    // Buffer parameters for tranforming 2D block-cyclic distributed DM matrix 
+// to grid distributed DM matrix
+    int* sender_2D_index;
+    int sender_size;
+    int* sender_size_process;
+    int* sender_displacement_process;
+    double* sender_buffer;
+
+    int* receiver_local_index;
+    int receiver_size;
+    int* receiver_size_process;
+    int* receiver_displacement_process;
+    double* receiver_buffer;
+#ifdef __MPI
+    MPI_Comm comm_2D;
+#endif
+};
\ No newline at end of file
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.cpp
index 2ad5f3dfa7e..aad613f17f0 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.cpp
@@ -30,15 +30,6 @@ Local_Orbital_Charge::Local_Orbital_Charge()
     // band_local = nullptr;
     // Z_wg = nullptr;
     // Z_LOC = nullptr;
-    sender_2D_index = nullptr;
-    sender_size_process = nullptr;
-    sender_displacement_process = nullptr;
-    sender_buffer = nullptr;
-
-    receiver_local_index = nullptr;
-    receiver_size_process = nullptr;
-    receiver_displacement_process = nullptr;
-    receiver_buffer = nullptr;
 }
 
 Local_Orbital_Charge::~Local_Orbital_Charge()
@@ -54,15 +45,6 @@ Local_Orbital_Charge::~Local_Orbital_Charge()
         delete[] DM;
         delete[] DM_pool;
     }
-    delete[] sender_2D_index;
-    delete[] sender_size_process;
-    delete[] sender_displacement_process;
-    delete[] sender_buffer;
-
-    delete[] receiver_local_index;
-    delete[] receiver_size_process;
-    delete[] receiver_displacement_process;
-    delete[] receiver_buffer;
 
     // with k points
     if (this->init_DM_R)
@@ -75,6 +57,7 @@ Local_Orbital_Charge::~Local_Orbital_Charge()
     }
 }
 
+
 void Local_Orbital_Charge::allocate_dm_wfc(const Grid_Technique &gt,
                                            elecstate::ElecState *pelec,
                                            Local_Orbital_wfc &lowf,
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h
index d8985bf6f8d..cfd58fce6b6 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h
@@ -12,7 +12,7 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h"
 #include "module_psi/psi.h"
 #include "module_elecstate/elecstate.h"
-
+#include "module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.h"
 class Local_Orbital_Charge
 {
 	public:
@@ -100,26 +100,10 @@ class Local_Orbital_Charge
 	// add by yshen on 9/22/2014
 	// these variables are memory pool for DM series matrixes, 
 	// so that these matrixes will be storaged continuously in the memory.
-	double **DM_pool;
-	
-	// Buffer parameters for tranforming 2D block-cyclic distributed DM matrix 
-	// to grid distributed DM matrix
-    int *sender_2D_index;
-    int sender_size;
-    int *sender_size_process;
-    int *sender_displacement_process;
-    double* sender_buffer;
-
-    int *receiver_local_index;
-    int receiver_size;
-    int *receiver_size_process;
-    int *receiver_displacement_process;
-    double* receiver_buffer;
-
-#ifdef __MPI
-    int setAlltoallvParameter(MPI_Comm comm_2D, int blacs_ctxt, int nblk);
-#endif
-    void cal_dk_gamma_from_2D(void);
+    double** DM_pool;
+
+    DMgamma_2dtoGrid dm2g;
+
 };
 
 #endif
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.cpp
index c9f1c8c0c12..248d01ddaaf 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.cpp
@@ -2,7 +2,6 @@
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/write_wfc_nao.h"
 #include "module_io/read_wfc_nao.h"
-#include "module_base/parallel_common.h"
 #include "module_base/memory.h"
 #include "module_base/timer.h"