diff --git a/source/module_cell/module_paw/paw_atom.cpp b/source/module_cell/module_paw/paw_atom.cpp
new file mode 100644
index 00000000000..61b285cfe81
--- /dev/null
+++ b/source/module_cell/module_paw/paw_atom.cpp
@@ -0,0 +1,60 @@
+#include "paw_atom.h"
+
+void Paw_Atom::init_paw_atom(const int nproj_in)
+{
+
+    nproj = nproj_in;
+
+    ca.resize(nproj);
+    rhoij.resize(nproj*(nproj + 1) / 2);
+    rhoijp.resize(nproj*(nproj + 1) / 2);
+    rhoijselect.resize(nproj*(nproj + 1) / 2);
+}
+
+void Paw_Atom::set_ca(std::vector<std::complex<double>> & ca_in, const double weight_in)
+{
+    for(int i = 0; i < nproj; i ++)
+    {
+        ca[i] = ca_in[i];
+    }
+
+    weight = weight_in;
+}
+
+void Paw_Atom::reset_rhoij()
+{
+    nrhoijsel = 0;
+    for(int i = 0; i < nproj*(nproj+1)/2; i ++)
+    {
+        rhoij[i] = 0.0;
+        rhoijp[i] = 0.0;
+        rhoijselect[i] = -1;
+    }    
+}
+
+void Paw_Atom::accumulate_rhoij()
+{
+    for(int iproj = 0; iproj < nproj; iproj ++)
+    {
+        int i0 = iproj * (iproj + 1) / 2;
+        for(int jproj = 0; jproj < iproj+1; jproj ++)
+        {
+            std::complex<double> tmp = std::conj(ca[iproj]) * ca[jproj];
+            rhoij[i0 + jproj] += tmp.real() * weight;
+        }
+    }
+}
+
+void Paw_Atom::convert_rhoij()
+{
+    nrhoijsel = 0;
+    for(int i = 0; i < rhoij.size(); i ++)
+    {
+        if(std::abs(rhoij[i]) > 1e-10)
+        {
+            rhoijselect[nrhoijsel] = i;
+            rhoijp[nrhoijsel] = rhoij[i];
+            nrhoijsel ++;
+        }
+    }
+}
\ No newline at end of file
diff --git a/source/module_cell/module_paw/paw_atom.h b/source/module_cell/module_paw/paw_atom.h
index fcc5ec22741..e8b2ec93d2b 100644
--- a/source/module_cell/module_paw/paw_atom.h
+++ b/source/module_cell/module_paw/paw_atom.h
@@ -5,27 +5,46 @@
 #define PAW_ATOM
 
 #include <vector>
+#include <complex>
 
 class Paw_Atom
 {
     public:
 
-    Paw_Atom();
-    ~Paw_Atom();
+    Paw_Atom(){};
+    ~Paw_Atom(){};
 
-    void calculate_rhoij(); //calculate <psi|ptilde><ptilde|psi>
+    //set the sizes of data structures
+    void init_paw_atom(const int nproj_in);
+
+    //pass <psi|ptilde> from outside and saves it
+    void set_ca(std::vector<std::complex<double>> & ca_in, const double weight_in);
+
+    void reset_rhoij(); //set rhoij = 0
+    void accumulate_rhoij(); //calculate and accumulate <psi|ptilde><ptilde|psi> from <psi|ptilde>
+
+    // not sure this is gonna be used, but it is nice to have
+    // an interface that returns rhoij and rhoijp I suppose
+    std::vector<double> get_rhoij(){return rhoij;}
+    std::vector<double> get_rhoijp(){return rhoijp;}
+    std::vector<int> get_rhoijselect(){return rhoijselect;}
+    int get_nrhoijsel(){return nrhoijsel;}
 
     void convert_rhoij(); //convert to format in libpaw
 
     private:
 
-    std::vector<std::vector<double>> ca(:,:) //coefficients <psi|ptilde>
-    std::vector<std::vector<double>> rhoij(:,:) //on-site density matrix
+    int nproj;
+
+    std::vector<std::complex<double>> ca; //coefficients <psi|ptilde> for a given psi
+    std::vector<double> rhoij; //on-site density matrix, upper triangular
+
+    double weight; //weight of current band
     
     //for libpaw
-    integer :: nrhoijsel; //#. nonzero elements
+    int nrhoijsel; //#. nonzero elements
     std::vector<int>    rhoijselect; //index of nonzero elements
-    std::vector<double> rhoij_paw;   //stores nonzero elements
+    std::vector<double> rhoijp; //rhoij in 'packed' format, only non-zero elements stored;
 };
 
 #endif
\ No newline at end of file
diff --git a/source/module_cell/module_paw/test/CMakeLists.txt b/source/module_cell/module_paw/test/CMakeLists.txt
index d8bda3c7e57..e3a5e6d1688 100644
--- a/source/module_cell/module_paw/test/CMakeLists.txt
+++ b/source/module_cell/module_paw/test/CMakeLists.txt
@@ -3,6 +3,7 @@ remove_definitions(-D__MPI)
 install(FILES sphbes_ref.dat func.dat qlist.dat fq_ref.dat
   qlist1.dat fq.dat d2fq_ref.dat qlist2.dat fq_fit_ref.dat
   kpg.dat ylm_ref.dat leg_pol.dat kpg1.dat igxyz.dat
+  ca.dat rhoij.dat
   H.LDA_PW-JTH.xml Fe.GGA_PBE-JTH.xml O.GGA_PBE-JTH.xml Si.GGA_PBE-JTH.xml
   DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 
@@ -16,4 +17,10 @@ AddTest(
   TARGET Test_Paw1
   LIBS ${math_libs} base device
   SOURCES test_paw1.cpp ../paw_element.cpp ../paw_sphbes.cpp ../paw_interface.cpp ../paw_cell.cpp
+)
+
+AddTest(
+  TARGET Test_Paw2
+  LIBS ${math_libs} base device
+  SOURCES test_paw2.cpp ../paw_atom.cpp
 )
\ No newline at end of file
diff --git a/source/module_cell/module_paw/test/ca.dat b/source/module_cell/module_paw/test/ca.dat
new file mode 100644
index 00000000000..06aa52e1d9f
--- /dev/null
+++ b/source/module_cell/module_paw/test/ca.dat
@@ -0,0 +1,24 @@
+  0.968906870089821       2.037955226635785E-002  7.396273824380536E-003
+ -2.042811836208747E-004  0.000000000000000E+000  0.000000000000000E+000
+ -9.289045159617711E-005  6.311929742912709E-003 -9.794059268907382E-005
+  2.365532747441629E-002  0.000000000000000E+000  0.000000000000000E+000
+ -1.142650727399384E-005  7.698728632265050E-004 -1.170958351262743E-005
+  2.828184157100050E-003
+ -7.033033646581735E-004 -9.303639153266802E-005  1.509561656125960E-004
+  9.449861517427804E-005  0.000000000000000E+000  0.000000000000000E+000
+  1.544265329194396E-003  8.419422128085135E-004  1.545133055206470E-002
+  4.949843069619822E-003  0.000000000000000E+000  0.000000000000000E+000
+  1.846259588735391E-004  1.020000215682333E-004  1.847330514478268E-003
+  5.917934454625722E-004
+  6.451295378324393E-004 -3.170596688672782E-005 -5.433598139695228E-005
+  2.139889904474025E-006  0.000000000000000E+000  0.000000000000000E+000
+ -2.417183344925036E-003 -3.175565864876030E-004 -2.390248730291169E-002
+ -4.202500707012095E-003  0.000000000000000E+000  0.000000000000000E+000
+ -2.890976125020449E-004 -3.757009410824718E-005 -2.857734097255318E-003
+ -5.024426710062551E-004
+  1.814549181115905E-003 -2.286152182260669E-004  1.208406080687989E-003
+ -7.097241086786346E-004  0.000000000000000E+000  0.000000000000000E+000
+ -4.641252446518110E-004  5.317470384156478E-004 -5.258243920811105E-003
+  6.174987459844208E-004  0.000000000000000E+000  0.000000000000000E+000
+ -5.525176964379127E-005  6.538856737945766E-005 -6.286652411423689E-004
+  7.382692851371659E-005
diff --git a/source/module_cell/module_paw/test/rhoij.dat b/source/module_cell/module_paw/test/rhoij.dat
new file mode 100644
index 00000000000..a2a058b260c
--- /dev/null
+++ b/source/module_cell/module_paw/test/rhoij.dat
@@ -0,0 +1,24 @@
+  1.87840022880052       1.432868449412428E-002  1.134904511673822E-004
+  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
+  6.990923125856316E-005 -4.942715832259250E-006  0.000000000000000E+000
+  9.876898495263231E-005  7.017854263446076E-004 -1.651816453291899E-005
+  0.000000000000000E+000  4.784567231285155E-004  2.879698035629865E-003
+  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
+  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
+  8.357169991407489E-006 -6.036424785199192E-007  0.000000000000000E+000
+  1.200516182219030E-005  5.793851926095698E-005  0.000000000000000E+000
+  1.459285808018694E-006  8.390407727891744E-005 -1.974879074783394E-006
+  0.000000000000000E+000  5.720333931853391E-005  3.442910004272299E-004
+  0.000000000000000E+000  6.927014746969510E-006  4.116275092338139E-005
+  1.87840022880053       1.432868449412449E-002  1.134904511673853E-004
+  6.990923125856056E-005 -4.942715832259143E-006  9.876898495263226E-005
+  7.017854263445924E-004 -1.651816453291776E-005  4.784567231285154E-004
+  2.879698035629866E-003  8.357169991407169E-006 -6.036424785199067E-007
+  1.200516182219029E-005  5.793851926095698E-005  1.459285808018692E-006
+  8.390407727891564E-005 -1.974879074783246E-006  5.720333931853392E-005
+  3.442910004272300E-004  6.927014746969510E-006  4.116275092338141E-005
+  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
+  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
+  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
+  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
+  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
diff --git a/source/module_cell/module_paw/test/test_paw2.cpp b/source/module_cell/module_paw/test/test_paw2.cpp
new file mode 100644
index 00000000000..4491484311c
--- /dev/null
+++ b/source/module_cell/module_paw/test/test_paw2.cpp
@@ -0,0 +1,69 @@
+#include "gtest/gtest.h"
+#include <fstream>
+#include <iostream>
+
+#include "../paw_atom.h"
+
+class Test_Paw_Atom : public testing::Test
+{
+    protected:
+
+    Paw_Atom paw_atom;
+    const int nproj = 8;
+};
+
+TEST_F(Test_Paw_Atom, test_paw)
+{
+    paw_atom.init_paw_atom(nproj);
+    paw_atom.reset_rhoij();
+
+    std::vector<std::complex<double>> ca;
+    ca.resize(nproj);
+    std::ifstream ifs_ca("ca.dat");
+    std::ifstream ifs_rhoij("rhoij.dat");
+
+    // there are altogether 4 bands
+    const int nband = 4;
+    for(int iband = 0; iband < nband; iband ++)
+    {
+        for(int ip = 0; ip < nproj; ip ++)
+        {
+            double re,im;
+            ifs_ca >> re >> im;
+            ca[ip] = std::complex<double>(re,im);
+        }
+
+        paw_atom.set_ca(ca, 2.0); //pass coefficient into the class
+        paw_atom.accumulate_rhoij(); //accumulate the contribution of current band
+    }
+
+    std::vector<double> rhoij = paw_atom.get_rhoij();
+
+    for(int i=0; i<rhoij.size();i++)
+    {
+        double tmp;
+        ifs_rhoij >> tmp;
+        EXPECT_NEAR(tmp,rhoij[i],1e-8);        
+    }
+
+    paw_atom.convert_rhoij();
+
+    int nrhoijsel = paw_atom.get_nrhoijsel();
+    EXPECT_EQ(nrhoijsel,21);
+
+    std::vector<int> rhoijselect = paw_atom.get_rhoijselect();
+    std::vector<int> rhoijselect_ref = {1,2,3,7,8,10,11,12,14,15,22,23,25,26,
+        28,29,30,32,33,35,36,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+    for(int i = 0; i < rhoijselect.size(); i ++)
+    {
+        EXPECT_EQ(rhoijselect[i],rhoijselect_ref[i]-1); // -1 because Fortran index starts from 1
+    }
+
+    std::vector<double> rhoijp = paw_atom.get_rhoijp();
+    for(int i=0; i<rhoijp.size();i++)
+    {
+        double tmp;
+        ifs_rhoij >> tmp;
+        EXPECT_NEAR(tmp,rhoijp[i],1e-8);        
+    }
+}
\ No newline at end of file