diff --git a/src/Decomposition/ie_ghost.hpp b/src/Decomposition/ie_ghost.hpp
index 8b1bbe5fbe219033d624b7d9edc0293f26850bb5..e6d65ffc4479dbe391a4592c1260f01d380a701b 100755
--- a/src/Decomposition/ie_ghost.hpp
+++ b/src/Decomposition/ie_ghost.hpp
@@ -31,11 +31,10 @@ class ie_ghost
 	//! It store the same information of box_nn_processor_int organized by processor id
 	openfpm::vector< Box_dom<dim,T> > proc_int_box;
 
-	// External ghost boxes for this processor, indicated with G8_0 G9_0 ...
+	//! External ghost boxes for this processor
 	openfpm::vector<p_box<dim,T> > vb_ext;
 
-	// Internal ghost boxes for this processor domain, indicated with B8_0 B9_0 ..... in the figure
-	// below as a linear vector
+	//! Internal ghost boxes for this processor domain
 	openfpm::vector<p_box<dim,T> > vb_int;
 
 	//! Cell-list that store the geometrical information of the internal ghost boxes
@@ -44,8 +43,10 @@ class ie_ghost
 	//! shift vectors
 	openfpm::vector<Point<dim,T>> shifts;
 
-	// Temporal buffers to return information for ghost_processorID
+	//! Temporal buffers to return temporal information for ghost_processorID
 	openfpm::vector<std::pair<size_t,size_t>> ids_p;
+
+	//! Temporal buffers to return temporal information
 	openfpm::vector<size_t> ids;
 
 
@@ -188,6 +189,9 @@ protected:
 	 * The geo cell list structure exist to speed up the labelling the points if they fall on some
 	 * internal ghost
 	 *
+	 * \param domain where the cell list is defined
+	 * \param div number of division of the cell list
+	 *
 	 */
 	void Initialize_geo_cell(const Box<dim,T> & domain, const size_t (&div)[dim])
 	{
@@ -200,7 +204,11 @@ protected:
 	 * For each sub-domain of the local processor it store the intersection between the enlarged
 	 * sub-domain of the calling processor with the adjacent processors sub-domains (External ghost box)
 	 *
+	 * \param v_cl Virtual cluster
 	 * \param ghost margins
+	 * \param subdomains vector of local sundomains
+	 * \param box_nn_processor it will store for each sub-domain the near processors
+	 * \param nn_prcs contain the sub-domains of the near processors
 	 *
 	 * \note Are the G8_0 G9_0 G9_1 G5_0 boxes in calculateGhostBoxes
 	 * \see calculateGhostBoxes
@@ -296,8 +304,8 @@ protected:
 	 * \param v_cl Virtual cluster
 	 * \param ghost margins
 	 * \param sub_domains
-	 * \param box_nn_processors sub-domains of the adjacent processors
-	 * \param nn_prcs structure that store the adjacent processor information
+	 * \param box_nn_processors sub-domains of the near processors
+	 * \param nn_prcs structure that store the near processor sub-domains
 	 * \param geo_cell Cell list that store the subdomain information
 	 *
 	 * \note Are the B8_0 B9_0 B9_1 B5_0 boxes in calculateGhostBoxes
@@ -508,6 +516,8 @@ public:
 	 * This function return the set of shift vectors that determine such shift, for example
 	 * in the example above the shift at position 5 will be (0,-1.0)
 	 *
+	 * \return the shift vectors
+	 *
 	 */
 	const openfpm::vector<Point<dim,T>> & getShiftVectors()
 	{
@@ -667,6 +677,8 @@ public:
 	}
 
 	/*! \brief Given the internal ghost box id, it return the internal ghost box
+	 *
+	 * \param b_id internal ghost box id
 	 *
 	 * \return the internal ghost box
 	 *
@@ -679,6 +691,8 @@ public:
 	/*! \brief Given the internal ghost box id, it return the near processor at witch belong
 	 *         or the near processor that produced this internal ghost box
 	 *
+	 * \param internal ghost box id
+	 *
 	 * \return the processor id of the ghost box
 	 *
 	 */
@@ -698,6 +712,8 @@ public:
 	}
 
 	/*! \brief Given the external ghost box id, it return the external ghost box
+	 *
+	 * \param b_id external ghost box id
 	 *
 	 * \return the external ghost box
 	 *
@@ -710,6 +726,8 @@ public:
 	/*! \brief Given the external ghost box id, it return the near processor at witch belong
 	 *         or the near processor that produced this external ghost box
 	 *
+	 * \param b_id external ghost box id
+	 *
 	 * \return the processor id of the external ghost box
 	 *
 	 */
@@ -721,6 +739,7 @@ public:
 	/*! /brief Given a point it return the set of boxes in which the point fall
 	 *
 	 * \param p Point to check
+	 *
 	 * \return An iterator with the id's of the internal boxes in which the point fall
 	 *
 	 */
@@ -757,6 +776,8 @@ public:
 	 *
 	 * \param return the processor ids (not the rank, the id in the near processor list)
 	 *
+	 * \return a vector of pairs containinf the requested infromation
+	 *
 	 */
 	template <typename id1, typename id2> inline const openfpm::vector<std::pair<size_t,size_t>> ghost_processorID_pair(Point<dim,T> & p, const int opt = MULTIPLE)
 	{
@@ -806,6 +827,8 @@ public:
 	 *
 	 * \param return the processor ids
 	 *
+	 * \return a vector containing the requested information
+	 *
 	 */
 	template <typename id> inline const openfpm::vector<size_t> ghost_processorID(Point<dim,T> & p, const int opt = MULTIPLE)
 	{
@@ -842,10 +865,12 @@ public:
 	/*! \brief Given a position it return if the position belong to any neighborhood processor ghost
 	 * (Internal ghost)
 	 *
-	 * \tparam id type of if to get box_id processor_id lc_processor_id
+	 * \tparam id1 first index type to get box_id processor_id lc_processor_id
+	 * \tparam id2 second index type to get box_id processor_id lc_processor_id
+	 *
 	 * \param p Particle position
 	 *
-	 * \param return the processor ids
+	 * \return a vector of pair containing the requested information
 	 *
 	 */
 	template<typename id1, typename id2, typename Mem> inline const openfpm::vector<std::pair<size_t,size_t>> ghost_processorID_pair(const encapc<1,Point<dim,T>,Mem> & p, const int opt = MULTIPLE)
diff --git a/src/Vector/vector_dist_cell_list_tests.hpp b/src/Vector/vector_dist_cell_list_tests.hpp
index abbd0e5c799577e4b1e933260328a39dd9004524..3826bbb9b74c49fb69a3cad570434fe8c682c5de 100644
--- a/src/Vector/vector_dist_cell_list_tests.hpp
+++ b/src/Vector/vector_dist_cell_list_tests.hpp
@@ -332,171 +332,6 @@ BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_reorder_forces_test )
 	}
 }
 
-/*
-BOOST_AUTO_TEST_CASE( vector_dist_sym_cell_list_test )
-{
-	long int k = 4096*create_vcluster().getProcessingUnits();
-
-	long int big_step = k / 30;
-	big_step = (big_step == 0)?1:big_step;
-	long int small_step = 21;
-
-	if (create_vcluster().getProcessingUnits() > 48)
-		return;
-
-	print_test( "Testing vector symmetric cell-list k<=",k);
-
-	// 3D test
-	for ( ; k > 8*big_step ; k-= (k > 2*big_step)?big_step:small_step )
-	{
-		double r_cut = 0.1;
-
-		// domain
-		Box<3,double> box({0.0,0.0,0.0},{1.0,1.0,1.0});
-
-		// Boundary conditions
-		size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
-
-		// ghost, big enough to contain the interaction radius
-		Ghost<3,float> ghost(r_cut);
-
-		vector_dist<3,double, aggregate<double> > vd(k,box,bc,ghost);
-
-		{
-		auto it = vd.getDomainIterator();
-
-		while (it.isNext())
-		{
-			auto p = it.get();
-
-			vd.getPos(p)[0] = (double)rand()/RAND_MAX;
-			vd.getPos(p)[1] = (double)rand()/RAND_MAX;
-			vd.getPos(p)[2] = (double)rand()/RAND_MAX;
-
-			++it;
-		}
-		}
-
-		vd.map();
-		vd.ghost_get<>();
-
-		// Get the Cell list structure
-		auto NN = vd.getCellList(r_cut);
-
-		// Get an iterator over particles
-		auto it2 = vd.getDomainAndGhostIterator();
-
-		openfpm::vector<openfpm::vector<size_t>> idx;
-		idx.resize(vd.size_local_with_ghost());
-
-		/////////// SYMMETRIC CASE CELL-LIST ////////
-
-		// For each particle ...
-		while (it2.isNext())
-		{
-			// ... p
-			auto p = it2.get();
-
-			// Get the position of the particle p
-			Point<3,double> xp = vd.getPos(p);
-
-			// Get an iterator over the neighborhood of the particle p symmetric
-			auto NpSym = NN.template getNNIteratorSym<NO_CHECK>(NN.getCell(vd.getPos(p)),p.getKey());
-
-			// For each neighborhood of the particle p
-			while (NpSym.isNext())
-			{
-				// Neighborhood particle q
-				auto q = NpSym.get();
-
-				// if p == q skip this particle
-				if (q == p.getKey() )	{++NpSym; continue;};
-
-				// Get position of the particle q
-				Point<3,double> xq = vd.getPos(q);
-
-				// take the normalized direction
-				double rn = norm2(xp - xq);
-
-				// potential energy (using pow is slower)
-				vd.getProp<0>(p) += rn;
-				vd.getProp<0>(q) += rn;
-
-				idx.get(p.getKey()).add(q);
-				idx.get(q).add(p.getKey());
-
-
-				// Next neighborhood
-				++NpSym;
-			}
-
-			// Next Particle
-			++it2;
-		}
-
-		/////////////// NON SYMMETRIC CASE ////////////////////////
-
-		openfpm::vector<openfpm::vector<size_t>> idx2;
-		idx2.resize(vd.size_local());
-
-		auto it = vd.getDomainIterator();
-
-		// For each particle ...
-		while (it.isNext())
-		{
-			// ... p
-			auto p = it.get();
-
-			// Get the position of the particle p
-			Point<3,double> xp = vd.getPos(p);
-
-			// Get an iterator over the neighborhood of the particle p
-			auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
-
-			double Ep = 0.0;
-
-			// For each neighborhood of the particle p
-			while (Np.isNext())
-			{
-				// Neighborhood particle q
-				auto q = Np.get();
-
-				// if p == q skip this particle
-				if (q == p.getKey())	{++Np; continue;};
-
-				// Get position of the particle q
-				Point<3,double> xq = vd.getPos(q);
-
-				// take the normalized direction
-				double rn = norm2(xp - xq);
-
-				idx2.get(p.getKey()).add(q);
-
-				// potential energy (using pow is slower)
-				Ep += rn;
-
-				// Next neighborhood
-				++Np;
-			}
-
-			idx.get(p.getKey()).sort();
-			idx2.get(p.getKey()).sort();
-
-			bool ret = true;
-
-			for (size_t i = 0 ; i < idx.get(p.getKey()).size() ; i++)
-				ret &= idx.get(p.getKey()).get(i) == idx2.get(p.getKey()).get(i);
-
-			BOOST_REQUIRE_EQUAL(ret,true);
-
-			BOOST_REQUIRE_CLOSE(Ep,vd.getProp<0>(p),0.01);
-
-			// Next Particle
-			++it;
-		}
-	}
-}*/
-
 BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 {
 	Vcluster & v_cl = create_vcluster();
@@ -504,11 +339,13 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 	if (v_cl.getProcessingUnits() > 24)
 		return;
 
+	float L = 1000.0;
+
     // set the seed
 	// create the random generator engine
-	std::srand(v_cl.getProcessUnitID());
+	std::srand(0);
     std::default_random_engine eg;
-    std::uniform_real_distribution<float> ud(0.0f, 1.0f);
+    std::uniform_real_distribution<float> ud(-L,L);
 
     long int k = 4096 * v_cl.getProcessingUnits();
 
@@ -518,20 +355,33 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 	print_test("Testing 3D periodic vector symmetric cell-list k=",k);
 	BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
 
-	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+	Box<3,float> box({-L,-L,-L},{L,L,L});
 
 	// Boundary conditions
 	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
-	float r_cut = 0.1;
+	float r_cut = 100.0;
 
 	// ghost
 	Ghost<3,float> ghost(r_cut);
 
-	typedef  aggregate<size_t,size_t> part_prop;
+	// Point and global id
+	struct point_and_gid
+	{
+		size_t id;
+		Point<3,float> xq;
+
+		bool operator<(const struct point_and_gid & pag)
+		{
+			return (id < pag.id);
+		}
+	};
+
+	typedef  aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
 
 	// Distributed vector
 	vector_dist<3,float, part_prop > vd(k,box,bc,ghost);
+	size_t start = vd.init_size_accum(k);
 
 	auto it = vd.getIterator();
 
@@ -545,7 +395,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 
 		// Fill some properties randomly
 
-		vd.getProp<0>(key) = 0.0;
+		vd.getProp<0>(key) = 0;
+		vd.getProp<1>(key) = 0;
+		vd.getProp<2>(key) = key.getKey() + start;
 
 		++it;
 	}
@@ -553,10 +405,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 	vd.map();
 
 	// sync the ghost
-	vd.ghost_get<0>();
-
-	auto NN = vd.getCellList(0.1);
+	vd.ghost_get<0,2>();
 
+	auto NN = vd.getCellList(r_cut);
 	auto p_it = vd.getDomainIterator();
 
 	while (p_it.isNext())
@@ -571,6 +422,12 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 		{
 			auto q = Np.get();
 
+			if (p.getKey() == q)
+			{
+				++Np;
+				continue;
+			}
+
 			// repulsive
 
 			Point<3,float> xq = vd.getPos(q);
@@ -581,7 +438,12 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 			// Particle should be inside 2 * r_cut range
 
 			if (distance < r_cut )
+			{
 				vd.getProp<0>(p)++;
+				vd.getProp<3>(p).add();
+				vd.getProp<3>(p).last().xq = xq;
+				vd.getProp<3>(p).last().id = vd.getProp<2>(q);
+			}
 
 			++Np;
 		}
@@ -591,7 +453,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 
 	// We now try symmetric  Cell-list
 
-	auto NN2 = vd.getCellListSym(0.1);
+	auto NN2 = vd.getCellListSym(r_cut);
 
 	auto p_it2 = vd.getDomainIterator();
 
@@ -607,6 +469,12 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 		{
 			auto q = Np.get();
 
+			if (p.getKey() == q)
+			{
+				++Np;
+				continue;
+			}
+
 			// repulsive
 
 			Point<3,float> xq = vd.getPos(q);
@@ -620,6 +488,14 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 			{
 				vd.getProp<1>(p)++;
 				vd.getProp<1>(q)++;
+
+				vd.getProp<4>(p).add();
+				vd.getProp<4>(q).add();
+
+				vd.getProp<4>(p).last().xq = xq;
+				vd.getProp<4>(q).last().xq = xp;
+				vd.getProp<4>(p).last().id = vd.getProp<2>(q);
+				vd.getProp<4>(q).last().id = vd.getProp<2>(p);
 			}
 
 			++Np;
@@ -629,6 +505,29 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 	}
 
 	vd.ghost_put<add_,1>();
+	vd.ghost_put<merge_,4>();
+
+	auto p_it3 = vd.getDomainIterator();
+
+	bool ret = true;
+	while (p_it3.isNext())
+	{
+		auto p = p_it3.get();
+
+		ret &= vd.getProp<1>(p) == vd.getProp<0>(p);
+
+		vd.getProp<3>(p).sort();
+		vd.getProp<4>(p).sort();
+
+		ret &= vd.getProp<3>(p).size() == vd.getProp<4>(p).size();
+
+		for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++)
+			ret &= vd.getProp<3>(p).get(i).id == vd.getProp<4>(p).get(i).id;
+
+		++p_it3;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
 }
 
 BOOST_AUTO_TEST_CASE( vector_dist_sym_verlet_list_test )
diff --git a/src/Vector/vector_dist_comm.hpp b/src/Vector/vector_dist_comm.hpp
index 0aae55388ac4b9aa559136d75f5f5b95e545222f..6a42f4c4658a0dd4ddd987c1d74744e54a1764be 100644
--- a/src/Vector/vector_dist_comm.hpp
+++ b/src/Vector/vector_dist_comm.hpp
@@ -1047,6 +1047,27 @@ public:
 		// Send and receive ghost particle information
 		op_ssend_recv_merge<op> opm(g_opart);
 		v_cl.SSendRecvP_op<op_ssend_recv_merge<op>,send_vector,decltype(v_prp),prp...>(g_send_prp,v_prp,prc_recv_get,opm,prc_recv_put,recv_sz_put);
+
+		// process also the local replicated particles
+
+		size_t i2 = 0;
+
+		#ifdef SE_CLASS1
+
+		if (v_prp.size() - lg_m != o_part_loc.size())
+			std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " v_prp.size() - lg_m = " << v_prp.size() - lg_m << " != " << o_part_loc.size() << std::endl;
+
+		#endif
+
+		for (size_t i = lg_m ; i < v_prp.size() ; i++)
+		{
+			auto dst = v_prp.get(o_part_loc.template get<0>(i2));
+			auto src = v_prp.get(i);
+			copy_cpu_encap_encap_op_prp<op,decltype(v_prp.get(0)),decltype(v_prp.get(0)),prp...> cp(src,dst);
+
+			boost::mpl::for_each_ref< boost::mpl::range_c<int,0,sizeof...(prp)> >(cp);
+			i2++;
+		}
 	}
 };
 
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index e76667dd2c3acc117614d672e5e70c75e1325eac..90bd5662832f41d3124dc50f7982fc543fcfd3fd 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -1618,13 +1618,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 	// 3D test
 	for ( ; k >= 2 ; k-= (k > 2*big_step)?big_step:small_step )
 	{
+		float r_cut = 1.3 / k;
+		float r_g = 1.5 / k;
+
 		Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
 
 		// Boundary conditions
 		size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
 		// ghost
-		Ghost<3,float> ghost(1.3/(k));
+		Ghost<3,float> ghost(r_g);
 
 		typedef  aggregate<float> part_prop;
 
@@ -1656,7 +1659,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 		vd.ghost_get<0>();
 
 		{
-			auto NN = vd.getCellList(1.3/k);
+			auto NN = vd.getCellList(r_cut);
 			float a = 1.0f*k*k;
 
 			// run trough all the particles + ghost
@@ -1680,8 +1683,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 
 					float dist = xp.distance(xq);
 
-					if (dist < 1.0/k)
-						vd.getProp<0>(q) += a*(-dist*dist+1.0/k/k);
+					if (dist < r_cut)
+						vd.getProp<0>(q) += a*(-dist*dist+r_cut*r_cut);
 
 					++Np;
 				}
@@ -1712,8 +1715,18 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 			BOOST_REQUIRE_EQUAL(ret,true);
 		}
 
+		auto itp = vd.getDomainAndGhostIterator();
+		while (itp.isNext())
+		{
+			auto key = itp.get();
+
+			vd.getProp<0>(key) = 0.0;
+
+			++itp;
+		}
+
 		{
-			auto NN = vd.getCellList(1.3/k);
+			auto NN = vd.getCellList(r_cut);
 			float a = 1.0f*k*k;
 
 			// run trough all the particles + ghost
@@ -1737,8 +1750,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 
 					float dist = xp.distance(xq);
 
-					if (dist < 1.0/k)
-						vd.getProp<0>(q) += a*(-dist*dist+1.0/k/k);
+					if (dist < r_cut)
+						vd.getProp<0>(q) += a*(-dist*dist+r_cut*r_cut);
 
 					++Np;
 				}