diff --git a/Src/Kernels/Chebyshev/FChebInterpolator.hpp b/Src/Kernels/Chebyshev/FChebInterpolator.hpp
index 97d525470235e9b73d98f95624776d97293da15c..e88439a48ef338bdbbc824cbc63ff3c4b695db59 100755
--- a/Src/Kernels/Chebyshev/FChebInterpolator.hpp
+++ b/Src/Kernels/Chebyshev/FChebInterpolator.hpp
@@ -637,14 +637,12 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
                                                                  FReal *const multipoleExpansion,
                                                                  const ContainerClass *const inParticles) const
 {
-    // set all multipole expansions to zero
-    // Should not be set to zero here (TODO remove this line) FBlas::setzero(2*nRhs*nnodes, multipoleExpansion);
 
     // allocate stuff
     const map_glob_loc map(center, width);
     FPoint localPosition;
 
-//    FReal W1 = FReal(0.);
+    // Init W
     FReal W1[nRhs];
     FReal W2[nRhs][3][ ORDER-1];
     FReal W4[nRhs][3][(ORDER-1)*(ORDER-1)];
@@ -750,7 +748,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
         for (unsigned int j=0; j<ORDER; ++j) {
             for (unsigned int k=0; k<ORDER; ++k) {
                 const unsigned int idx = k*ORDER*ORDER + j*ORDER + i;
-                multipoleExpansion[2*idxRhs*nnodes + idx] = (W1[idxRhs] +
+                multipoleExpansion[2*idxRhs*nnodes + idx] += (W1[idxRhs] +
                                                                      FReal(2.) * (F2[0][i] + F2[1][j] + F2[2][k]) +
                                                                      FReal(4.) * (F4[0][j*ORDER+i] + F4[1][k*ORDER+i] + F4[2][k*ORDER+j]) +
                                                                      FReal(8.) *  F8[idx]) / nnodes; // 11 * ORDER*ORDER*ORDER flops
diff --git a/Src/Kernels/Uniform/FUnifInterpolator.hpp b/Src/Kernels/Uniform/FUnifInterpolator.hpp
index cac8737764377a567e24216153ea1f93099b5f7c..9f63e6207585a63cf16f00b69bec062e0dea4ba2 100755
--- a/Src/Kernels/Uniform/FUnifInterpolator.hpp
+++ b/Src/Kernels/Uniform/FUnifInterpolator.hpp
@@ -415,8 +415,6 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
                                                                  FReal *const multipoleExpansion,
                                                                  const ContainerClass *const inParticles) const
 {
-  // set all multipole expansions to zero
-  // Should not be set to zero here (TODO remove this line) FBlas::setzero(nRhs*nnodes, multipoleExpansion);
 
   // allocate stuff
   const map_glob_loc map(center, width);