diff --git a/CHANGELOG.md b/CHANGELOG.md
index bfc246e8c1fc3922df0c97957e3d071df108df0b..7c85736920c4e3e09a0ac5cbba09e11771300269 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,12 +1,30 @@
 # Change Log
 All notable changes to this project will be documented in this file.
 
+## [0.6.0] - End October 2016
+
+### Added
+- Symmetric cell-list/verlet list
+- Multi-phase cell-list and Multi-phase cell-list
+- Added ghost_get that keep properties
+- Examples: 1_ghost_get_put it show how to use ghost_get and put with the new options
+            4_multiphase_celllist_verlet completely rewritten for new Cell-list and multiphase verlet
+	    5_molecular_dynamic use case of symmetric cell-list and verlet list with ghost put
+	    6_complex_usage It show how the flexibility of openfpm can be used to debug your program
+
+### Fixed
+- Option NO_POSITION was untested
+
+### Changes
+
+
 ## [0.5.1] - 27 September 2016
 
 ### Added
 - ghost_put support for particles
 - Full-Support for complex property on vector_dist (Serialization)
 - Added examples for serialization of complex properties 4_Vector
+- improved speed of the iterators
 
 ### Fixed
 - Installation PETSC installation fail in case of preinstalled MPI
@@ -131,17 +149,14 @@ All notable changes to this project will be documented in this file.
 
 - Algebraic Multigrid solver
 - Parallel VTK, improved visualization
+- Asynchronous communication
 
-## [0.6.0] - Middle of October
+## [0.7.0] - December of October
 
 ### Added
+- Asynchronous communication
 
-- Symmetric Cell list and Verlet (15 days)
-- Semantic communication (??)
-- Improved Finite difference interface (15 days) 
 
 
-## [0.6.0] - Beginning of September
 
-- Complex properties and serialization interface (15 days)
 
diff --git a/example/Vector/0_simple/main.cpp b/example/Vector/0_simple/main.cpp
index aeba11f0ade5818a5f0525908848fe4b26d8e959..0e6a334477539f0632710a74a26b19a9236587e4 100644
--- a/example/Vector/0_simple/main.cpp
+++ b/example/Vector/0_simple/main.cpp
@@ -348,7 +348,7 @@ int main(int argc, char* argv[])
 	 *
 	 * ## Finalize ## {#finalize_e0_sim}
 	 *
-	 *  At the very end of the program we have always to de-initialize the library
+	 *  At the very end of the program we have always de-initialize the library
 	 *
 	 * \snippet Vector/0_simple/main.cpp finalize
 	 *
diff --git a/example/Vector/4_multiphase_celllist_verlet/Makefile b/example/Vector/4_multiphase_celllist_verlet/Makefile
index 16168a458fb8b2d9c37bbe9e01d0b1ed7a3fbbaf..6043bb602cf6e5011ede7a7379252105c4d91db6 100644
--- a/example/Vector/4_multiphase_celllist_verlet/Makefile
+++ b/example/Vector/4_multiphase_celllist_verlet/Makefile
@@ -9,16 +9,16 @@ OBJ = main.o
 %.o: %.cpp
 	$(CC) -fext-numeric-literals -O3 -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
 
-cell: $(OBJ)
+multip: $(OBJ)
 	$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
 
-all: cell
+all: multip
 
 run: all
-	source $$HOME/openfpm_vars; mpirun -np 2 ./cell
+	source $$HOME/openfpm_vars; mpirun -np 2 ./multip
 
 .PHONY: clean all run
 
 clean:
-	rm -f *.o *~ core cell
+	rm -f *.o *~ core multip
 
diff --git a/example/Vector/4_multiphase_celllist_verlet/main.cpp b/example/Vector/4_multiphase_celllist_verlet/main.cpp
index c353f0d8d9897554c23b7059eae67e5350e17ff5..bf6223ca414c8bd0306cb574606811c2f0b9bb5a 100644
--- a/example/Vector/4_multiphase_celllist_verlet/main.cpp
+++ b/example/Vector/4_multiphase_celllist_verlet/main.cpp
@@ -3,16 +3,18 @@
 #include "Decomposition/CartDecomposition.hpp"
 #include "data_type/aggregate.hpp"
 #include "NN/CellList/CellListM.hpp"
+#include "Vector/vector_dist_multiphase_functions.hpp"
 
 /*!
- * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+ * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list and verlet
  *
  * [TOC]
  *
  *
- * # Vector Multi Phase cell-list # {#e4_ph_cl}
+ * # Vector Multi Phase cell-list and Verlet # {#e4_ph_cl}
  *
- * This example show multi-phase cell lists for the distributed vector
+ * This example show how to use multi-phase cell lists and Verlet-list.More in general
+ * it show how to construct Verlet and Cell-list between multiple vector_dist.
  *
  *
  */
@@ -25,19 +27,23 @@ int main(int argc, char* argv[])
 	 * ## Initialization ##
 	 *
 	 * Here we Initialize the library, and we create a set of distributed vectors all forced to have the same
-	 * decomposition. Each vector identify one phase
+	 * decomposition. Each vector identify one phase.
 	 *
-	 * \snippet Vector/4_multiphase_celllist/main.cpp Initialization and parameters
+	 * \note Be carefull on how you initialize the other phases. All the other phases
+	 *       must be forced to use the same decomposition. In order to do this we have
+	 *       to use the special constructor where we pass the decomposition from the first
+	 *       phase. The second parameter is just the the number of particles
+	 *
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp Initialization and parameters
 	 *
 	 */
 
 	//! \cond [Initialization and parameters] \endcond
 
 	openfpm_init(&argc,&argv);
-	Vcluster & v_cl = create_vcluster();
 
-	// we will place the particles on a grid like way with 128 particles on each direction
-	size_t sz[3] = {128,128,128};
+	// Vcluster for general usage
+	Vcluster & v_cl = create_vcluster();
 
 	// The domain
 	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
@@ -45,11 +51,13 @@ int main(int argc, char* argv[])
 	// Boundary conditions
 	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
+	// cut-off radius
 	float r_cut = 0.05;
 
 	// ghost, big enough to contain the interaction radius
 	Ghost<3,float> ghost(r_cut);
 
+	// The set of phases
 	openfpm::vector< vector_dist<3,float, aggregate<double,double>> > phases;
 
 	// first phase
@@ -67,108 +75,431 @@ int main(int argc, char* argv[])
 	/*!
 	 * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
 	 *
-	 * ## Initialization ##
+	 * ## Create phases ##
+	 *
+	 * We initialize all the phases with particle randomly positioned in the space. Completed this
+	 * iteration we redistribute the particles using the classical map, and we synchronize the ghost
 	 *
-	 * We initialize all the phases with particle randomly positioned in the space
+	 * \warning we cannot use the same iterator for all the phases even if the number of particles for
+	 *          each phase is the same. Consider that the number of particles (per-processor) can be
+	 *          different. The case of the initialization ("before map") is the only case where we are
+	 *          sure that the number of particles per processor is the same for each phase
 	 *
-	 * \snippet Vector/4_multiphase_celllist/main.cpp rand dist
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp rand dist
 	 *
 	 */
 
 	//! \cond [rand dist] \endcond
 
+	// An iterator over the particles of the phase0
 	auto it = phases.get(0).getDomainIterator();
 
-	// For all the particles
+	// For all the particles of phase0
 	while (it.isNext())
 	{
-		// for all phases
-		for (size_t i = 0 ; i < phases.size() ; i++)
-		{
-			auto key = it.get();
+		// particle p
+		auto p = it.get();
 
-			phases.get(i).getPos(key)[0] = (float)rand() / RAND_MAX;
-			phases.get(i).getPos(key)[1] = (float)rand() / RAND_MAX;
-			phases.get(i).getPos(key)[2] = (float)rand() / RAND_MAX;
+		// Assign the position of the particles to each phase
+		for (size_t i = 0 ; i < phases.size(); i++)
+		{
+			// Assign random position
+			phases.get(i).getPos(p)[0] = (float)rand() / RAND_MAX;
+			phases.get(i).getPos(p)[1] = (float)rand() / RAND_MAX;
+			phases.get(i).getPos(p)[2] = (float)rand() / RAND_MAX;
 		}
-		// next point
+
+		// Next particle
 		++it;
 	}
 
-	// Redistribute all the phases in the mean while also take the iterators of all the phases
 
-	typedef decltype(phases.get(0).getIterator()) iterator;
-	openfpm::vector<iterator> phase_it;
-
-	for (size_t i = 0 ; i < phases.size() ; i++)
+	// Redistribute and sync all the phases
+	for (size_t i = 0 ; i < 4 ; i++)
+	{
 		phases.get(i).map();
+		phases.get(i).ghost_get<>();
+	}
 
 	//! \cond [rand dist] \endcond
 
 	/*!
 	 * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
 	 *
-	 * ## Multi-phase cell-list construction ##
+	 * ## Multi phase Verlet-list ##
+	 *
+	 * In general verlet-list can be constructed from the vector itself using **getVerlerList()**
+	 *
+	 * \see \ref Vector_3_md_vl
+	 *
+	 * In the multi-phase case if we use such function on phase0 such function
+	 * produce a Verlet list where for each particle of the phase0 we get the neighborhood
+	 *  particles within the phase0. Most of time we also want to construct Verlet-list across
+	 *  phases, for example between phase0 and phase1. Let suppose now that we want to construct
+	 *  a verlet-list that given the particles of the phase0 it return the neighborhood particles
+	 *   in phase1. In order to do this we can create a Cell-list from phase1. Once we have the
+	 *   Cell-list for phase 1 we give to **createVerlet()**
+	 *  the particle set of phase0, the Cell-list of phase1 and the cut-off radius.
+	 *  The function return a VerletList between phase0 and 1 that can be used to do computation.
+	 *
+	 *
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp create multi-phase verlet
+	 *
+	 */
+
+	//! \cond [create multi-phase verlet] \endcond
+
+	// Get the cell list of the phase1
+	auto CL_phase1 = phases.get(1).getCellList(r_cut);
+
+	// This function create a Verlet-list between phases 0 and 1
+	auto NN_ver01 = createVerlet(phases.get(0),CL_phase1,r_cut);
+
+	//! \cond [create multi-phase verlet] \endcond
+
+	/*!
+	 * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+	 *
+	 * Once we have a Verlet-list, we can do easily computation with that. We can use
+	 * the function **getNNIterator()** to get an iterator of the neighborhood particles
+	 *  for a specified particle. In this case for each particle of the phase0 we count
+	 *  the neighborhood particles in phase1
+	 *
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp count part from phase0 to 1
+	 *
+	 */
+
+	//! \cond [count part from phase0 to 1] \endcond
+
+	// Get an iterator of the particles of the phase0
+	auto it = phases.get(0).getDomainIterator();
+
+	// For each particle of the phase0
+	while (it.isNext())
+	{
+		// Get the particle p
+		auto p = it.get();
+
+		// reset the counter
+		phases.get(0).getProp<0>(p) = 0;
+
+		// Get an iterator over the neighborhood particles for the particle p
+		auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey());
+
+		// For each neighborhood of the particle p
+		while (Np.isNext())
+		{
+			// Neighborhood particle q
+			auto q = Np.get();
+
+			// Count the number of particles
+			// Here in general we can do our computation
+			phases.get(0).getProp<0>(p)++;
+
+			// Next particle
+			++Np;
+		}
+
+		// Next particle
+		++it;
+	}
+
+	//! \cond [count part from phase0 to 1] \endcond
+
+	/*!
+	 * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+	 *
+	 * ## From one phase to all ##
+	 *
+	 * It is also possible to construct a Verlet-list from one phase to all the other phases.
+	 * To do this we use the function **createCellListM<2>()** to first create a
+	 * multi-phase cell-list. The template parameter is required to indicate how many bit
+	 * to reserve for the phase information.
+	 *
+	 * \note In case of
+	 * 	* 2 bit mean that you can store up to 4 phases (2^62 number of particles for each phase)
+	 * 	* 3 bit mean that you can store up to 8 phases (2^61 number of particles for each phase)
 	 *
-	 * In this part we construct the Multi-phase cell list. The multiphase cell list has 3 parameters
-	 * * one is the dimensionality (3)
-	 * * The precision of the coordinates (float),
-	 * * How many bit to use for the phase.
-	 *   Multi-phase cell-list try to pack into a 64bit number information about the particle id and the
-	 *   the phase id.
+	 * Once created the multiphase-cell list from the phases. We can use the function
+	 * **createVerletM()** to create a verlet-list from one phase to all the others.
+	 * The function requires the particle for which we are constructing the verlet list, in this case
+	 * the phase0, the Cell-list containing all the other phases and the cut-off radius.
 	 *
-	 * \snippet Vector/4_multiphase_celllist/main.cpp cl construction
+	 *
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp create multi-phase multi verlet
 	 *
 	 */
 
-	//! \cond [cl construction] \endcond
+	//! \cond [create multi-phase multi verlet] \endcond
+
+	// This function create an "Empty" Multiphase Cell List from all the phases
+	auto CL_all = createCellListM<2>(phases,r_cut);
+
+	// This create a Verlet-list between phase0 and all the other phases
+	auto NNver0_all = createVerletM<2>(phases.get(0),CL_all,r_cut);
 
-	// Construct one single Multi-phase cell list to use in the computation
-	// in 3d, precision float, 2 bit dedicated to the phase for a maximum of 2^2 = 4 (Maximum number of phase)
-	//
-	//
-	
-	size_t div[3];
-	Box<3,float> box_cl;
-	phases.get(0).getCellListParams(r_cut,div,box_cl);
+	//! \cond [create multi-phase multi verlet] \endcond
 
-	CellListM<3,float,2> NN;
-	NN.Initialize(box_cl,div);
+	/*!
+	 * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+	 *
+	 * Compute on a multiphase verlet-list is very similar to a non multi-phase.
+	 * The only difference is that the function **get()** is substituted by two
+	 * other functions **getP()** and **getV()** The first return the phase from where
+	 * it come the particle the second it return the particle-id
+	 *
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp compute multi-phase multi verlet
+	 *
+	 */
 
+	//! \cond [compute multi-phase multi verlet] \endcond
 
+	// Get an iterator for the phase0 particles
+	it = phases.get(0).getDomainIterator();
 
-	// for all the phases i
-	for (size_t i = 0; i < phases.size() ; i++)
+	while (it.isNext())
 	{
-		// iterate across all the particle of the phase i
-		auto it = phases.get(i).getDomainIterator();
+		// Get the particle p
+		auto p = it.get();
+
+		// Get an interator over the neighborhood of the particle p
+		auto Np = NNver0_all.template getNNIterator<NO_CHECK>(p.getKey());
+
+		// reset the counter
+		phases.get(0).getProp<0>(p) = 0;
+
+		// For each neighborhood of the particle p
+		while (Np.isNext())
+		{
+			// Get the particle q near to p
+			auto q = Np.getP();
+
+			// Get from which phase it come from
+			auto ph_q = Np.getV();
+
+			// count
+			phases.get(0).getProp<0>(p)++;
+
+			// Next particle
+			++Np;
+		}
+
+		// Next particle
+		++it;
+	}
+
+	//! \cond [compute multi-phase multi verlet] \endcond
+
+	/*!
+	 * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+	 *
+	 * ## Symmetric interaction case ##
+	 *
+	 * The same functions exist also in the case we want to construct
+	 * symmetric-verlet list.
+	 *
+	 * \see \ref Vector_5_md_vl_sym For more details on how to use symmetric
+	 * verlet-list in a real case.
+	 *
+	 * In general the main differences can be summarized in
+	 * * we have to reset the **forces** or **interaction** variable(s)
+	 * * calculate the interaction p-q and apply the result to p and q at the same time
+	 * * merge with ghost_put the result is stored in the ghost part with the
+	 *   real particles
+	 *
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp compute sym multi-phase two phase
+	 *
+	 */
+
+	//! \cond [compute sym multi-phase two phase] \endcond
+
+	// Get the cell list of the phase1
+	auto CL_phase1 = phases.get(1).getCellListSym(r_cut);
+
+	// This function create a Verlet-list between phases 0 and 1
+	auto NN_ver01 = createVerletSym(phases.get(0),CL_phase1,r_cut);
+
+	// Get an iterator over the real and ghost particles
+	auto it = phases.get(0).getDomainAndGhostIterator();
+
+	// For each particles
+	while (it.isNext())
+	{
+		// Get the particle p
+		auto p = it.get();
+
+		// Reset the counter
+		phases.get(0).getProp<0>(p) = 0;
+
+		// Next particle
+		++it;
+	}
+
+	// Compute interaction from phase0 to phase1
+
+	// Get an iterator over the real particles of phase0
+	auto it = phases.get(0).getDomainIterator();
+
+	// For each particle of the phase0
+	while (it.isNext())
+	{
+		// get particle p
+		auto p = it.get();
+
+		// Get the neighborhood of particle p
+		auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey());
+
+		// For each neighborhood of the particle p
+		while (Np.isNext())
+		{
+			// Neighborhood particle q of p
+			auto q = Np.get();
+
+			// increment the counter for the phase 0 and 1
+			phases.get(0).getProp<0>(p)++;
+			phases.get(1).getProp<0>(q)++;
+
+			// Next particle
+			++Np;
+		}
+
+		// Next particle
+		++it;
+	}
+
+	// Merge the information of the ghost with the real particles
+	phases.get(0).ghost_put<add_,0>();
+	phases.get(1).ghost_put<add_,0>();
+
+	//! \cond [compute sym multi-phase two phase] \endcond
+
+	/*!
+	 * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
+	 *
+	 * For the case of a Verlet-list between the phase0 and all the other phases.
+	 * The code remain very similar to the previous one the only differences is that
+	 * we have to create a Multi-phase cell list from the vector of the phases.
+	 * When we compute with the verlet-list list now the simple function **get**
+	 * is substituted by the function **getP** and **getV**. In this example we are
+	 * creating 4 multiphase symmetric verlet-list
+	 *
+	 * * 0 to all
+	 * * 1 to all
+	 * * 2 to all
+	 * * 3 to all
+	 *
+	 * The computation is an all to all
+	 *
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp create sym multi-phase multi verlet
+	 *
+	 */
+
+	//! \cond [create sym multi-phase multi verlet] \endcond
+
+	// Get an iterator over the phase0
+	auto it = phases.get(0).getDomainAndGhostIterator();
+
+	// For each particle of the phase 0
+	while (it.isNext())
+	{
+		// get the particle p
+		auto p = it.get();
+
+		// Reset the counter for the domain and ghost particles
+		phases.get(0).getProp<0>(p) = 0;
+		phases.get(1).getProp<0>(p) = 0;
+		phases.get(2).getProp<0>(p) = 0;
+		phases.get(3).getProp<0>(p) = 0;
+
+		// next particle
+		++it;
+	}
+
+	// This function create an "Empty" Multiphase Cell List
+	auto CL_all = createCellListSymM<2>(phases,r_cut);
+
+	// Type of the multiphase Verlet-list
+	typedef decltype(createVerletSymM<2>(phases.get(0),CL_all,r_cut)) verlet_type;
+
+	// for each phase we create one Verlet-list that contain the neighborhood
+	// from all the phases
+	verlet_type NNver_all[4];
+
+	// Here we create a Verlet-list between each phases
+
+	// 0 to all
+	NNver_all[0] = createVerletSymM<2>(phases.get(0),CL_all,r_cut);
+
+	// 1 to all
+	NNver_all[1] = createVerletSymM<2>(phases.get(1),CL_all,r_cut);
+
+	// 2 to all
+	NNver_all[2] = createVerletSymM<2>(phases.get(2),CL_all,r_cut);
+
+	// 3 to all
+	NNver_all[3] = createVerletSymM<2>(phases.get(3),CL_all,r_cut);
+
+	// For each phase
+	for (size_t i = 0 ; i < phases.size() ; i++)
+	{
+		// Get an iterator over the particles of that phase
+		it = phases.get(i).getDomainIterator();
+
+		// for each particle of the phase
 		while (it.isNext())
 		{
-			auto key = it.get();
+			// Get the particle p
+			auto p = it.get();
+
+			// Get an iterator for neighborhood of the particle p
+			auto Np = NNver_all[i].template getNNIterator<NO_CHECK>(p.getKey());
+
+			// For each neighborhood particle
+			while (Np.isNext())
+			{
+				// Get the particle q near to p
+				auto q = Np.getP();
 
-			// Add the particle of the phase i to the cell list
-			NN.add(phases.get(i).getPos(key), key.getKey(), i);
+				// Get from which phase it come from
+				auto ph_q = Np.getV();
 
+				// increment the counter on both p and q
+				phases.get(i).getProp<0>(p)++;
+				phases.get(ph_q).getProp<0>(q)++;
+
+				// Next particle
+				++Np;
+			}
+
+			// Next particle
 			++it;
 		}
 	}
 
-	//! \cond [cl construction] \endcond
+	// Merge the ghost part with the real particles
+	phases.get(0).ghost_put<add_,0>();
+	phases.get(1).ghost_put<add_,0>();
+	phases.get(2).ghost_put<add_,0>();
+	phases.get(3).ghost_put<add_,0>();
+
+	//! \cond [create sym multi-phase multi verlet] \endcond
 
 	/*!
 	 * \page Vector_4_mp_cl Vector 4 Multi Phase cell-list
 	 *
 	 * ## Multi-phase cell-list usage ##
 	 *
-	 * After construction we show how to use the Cell-list. In this case we accumulate on the property
+	 * The multiphase cell-list that we constructed before to create a Verlet-list can be also used
+	 * directly. In this case we accumulate on the property
 	 * 0 of the phase 0 the distance of the near particles from all the phases
 	 *
-	 * \snippet Vector/4_multiphase_celllist/main.cpp cl usage
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp cl usage
 	 *
 	 */
 
 	//! \cond [cl usage] \endcond
 
+	// we create a reference to phase0 particle for convenience
 	vector_dist<3,float, aggregate<double,double> > & current_phase = phases.get(0);
 
 	// Get the iterator of the particles of phase 0
@@ -183,8 +514,11 @@ int main(int argc, char* argv[])
 		// Get the position of the particle p
 		Point<3,float> xp = current_phase.getPos(p);
 
+		// Reset to zero the propety 0
+		current_phase.getProp<0>(p) = 0.0;
+
 		// Get an iterator of all the particles neighborhood of p
-		auto Np = NN.getNNIterator(NN.getCell(current_phase.getPos(p)));
+		auto Np = CL_all.getNNIterator(NN.getCell(current_phase.getPos(p)));
 
 		// For each particle near p
 		while (Np.isNext())
@@ -216,7 +550,7 @@ int main(int argc, char* argv[])
 	 *
 	 *  At the very end of the program we have always de-initialize the library
 	 *
-	 * \snippet Vector/4_multiphase_celllist/main.cpp finalize
+	 * \snippet Vector/4_multiphase_celllist_verlet/main.cpp finalize
 	 *
 	 */
 
@@ -231,7 +565,7 @@ int main(int argc, char* argv[])
 	 *
 	 * # Full code # {#code}
 	 *
-	 * \include Vector/4_multiphase_celllist/main.cpp
+	 * \include Vector/4_multiphase_celllist_verlet/main.cpp
 	 *
 	 */
 }
diff --git a/src/Grid/staggered_dist_grid_util.hpp b/src/Grid/staggered_dist_grid_util.hpp
index c7d20ec6a1902d11264aa8ceeb7c67119ff368f8..938f1d9ab473645b22c4025043b49c17df9567fa 100644
--- a/src/Grid/staggered_dist_grid_util.hpp
+++ b/src/Grid/staggered_dist_grid_util.hpp
@@ -66,11 +66,13 @@ struct vtk_write<ele,vtk,false>
 template<typename T>
 struct extends
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return 1;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 0;
@@ -81,11 +83,13 @@ struct extends
 template<typename T,size_t N1>
 struct extends<T[N1]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 1;
@@ -96,11 +100,13 @@ struct extends<T[N1]>
 template<typename T,size_t N1,size_t N2>
 struct extends<T[N1][N2]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1 * N2;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 2;
@@ -111,11 +117,13 @@ struct extends<T[N1][N2]>
 template<typename T,size_t N1,size_t N2,size_t N3>
 struct extends<T[N1][N2][N3]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1 * N2 * N3;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 3;
@@ -126,11 +134,13 @@ struct extends<T[N1][N2][N3]>
 template<typename T,size_t N1,size_t N2,size_t N3,size_t N4>
 struct extends<T[N1][N2][N3][N4]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1 * N2 * N3 * N4;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 4;
@@ -141,11 +151,13 @@ struct extends<T[N1][N2][N3][N4]>
 template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5>
 struct extends<T[N1][N2][N3][N4][N5]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1 * N2 * N3 * N4 * N5;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 5;
@@ -156,11 +168,13 @@ struct extends<T[N1][N2][N3][N4][N5]>
 template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6>
 struct extends<T[N1][N2][N3][N4][N5][N6]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1 * N2 * N3 * N4 * N5 * N6;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 6;
@@ -171,11 +185,13 @@ struct extends<T[N1][N2][N3][N4][N5][N6]>
 template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7>
 struct extends<T[N1][N2][N3][N4][N5][N6][N7]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1 * N2 * N3 * N4 * N5 * N6 * N7;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 7;
@@ -186,11 +202,13 @@ struct extends<T[N1][N2][N3][N4][N5][N6][N7]>
 template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8>
 struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 8;
@@ -201,11 +219,13 @@ struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8]>
 template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9>
 struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 9;
@@ -216,11 +236,13 @@ struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9]>
 template<typename T,size_t N1,size_t N2,size_t N3,size_t N4,size_t N5, size_t N6, size_t N7, size_t N8, size_t N9, size_t N10>
 struct extends<T[N1][N2][N3][N4][N5][N6][N7][N8][N9][N10]>
 {
+	//! number of elements
 	static inline size_t mul()
 	{
 		return N1 * N2 * N3 * N4 * N5 * N6 * N7 * N8 * N9 * N10;
 	}
 
+	//! number of indexes
 	static inline size_t dim()
 	{
 		return 10;
diff --git a/src/Makefile.am b/src/Makefile.am
index fff4127b0a3b696e5dc7eed7d653e86ddf9b38b1..3ba89adcd433d3aaf5eca6a0b02af994b68e499f 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -9,7 +9,7 @@ nobase_include_HEADERS = Decomposition/CartDecomposition.hpp Decomposition/CartD
          Decomposition/nn_processor.hpp Decomposition/ie_loc_ghost.hpp Decomposition/ORB.hpp \
          Graph/CartesianGraphFactory.hpp \
          Grid/grid_dist_id.hpp Grid/grid_dist_id_iterator_dec.hpp Grid/grid_dist_util.hpp  Grid/grid_dist_id_iterator_sub.hpp Grid/grid_dist_id_iterator.hpp Grid/grid_dist_key.hpp Grid/staggered_dist_grid.hpp Grid/staggered_dist_grid_util.hpp Grid/staggered_dist_grid_copy.hpp \
-         Vector/vector_dist_comm.hpp Vector/vector_dist.hpp Vector/vector_dist_ofb.hpp Vector/vector_dist_iterator.hpp Vector/vector_dist_key.hpp \
+         Vector/vector_dist_multiphase_functions.hpp Vector/vector_dist_comm.hpp Vector/vector_dist.hpp Vector/vector_dist_ofb.hpp Vector/vector_dist_iterator.hpp Vector/vector_dist_key.hpp \
          config/config.h \
          example.mk \
          Decomposition/Distribution/metis_util.hpp Decomposition/Distribution/parmetis_dist_util.hpp  Decomposition/Distribution/parmetis_util.hpp Decomposition/Distribution/MetisDistribution.hpp Decomposition/Distribution/ParMetisDistribution.hpp Decomposition/Distribution/DistParMetisDistribution.hpp  dec_optimizer.hpp SubdomainGraphNodes.hpp \
diff --git a/src/Vector/vector_dist_MP_unit_tests.hpp b/src/Vector/vector_dist_MP_unit_tests.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7d9eff81e52f5cbb5b5e2ad3b1951e69802ffbe0
--- /dev/null
+++ b/src/Vector/vector_dist_MP_unit_tests.hpp
@@ -0,0 +1,377 @@
+/*
+ * vector_dist_MP_unit_tests.hpp
+ *
+ *  Created on: Oct 14, 2016
+ *      Author: i-bird
+ */
+
+#ifndef SRC_VECTOR_VECTOR_DIST_MP_UNIT_TESTS_HPP_
+#define SRC_VECTOR_VECTOR_DIST_MP_UNIT_TESTS_HPP_
+
+#include "Vector/vector_dist_multiphase_functions.hpp"
+
+BOOST_AUTO_TEST_SUITE( vector_dist_multiphase_test )
+
+BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_test )
+{
+	if (create_vcluster().getProcessingUnits() > 24)
+		return;
+
+	size_t sz[3] = {60,60,40};
+
+	// The domain
+	Box<3,float> box({-1000.0,-1000.0,-1000.0},{2000.0,2000.0,1000.0});
+
+	// Boundary conditions
+	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+	float r_cut = 51.0;
+
+	// ghost, big enough to contain the interaction radius
+	Ghost<3,float> ghost(r_cut);
+
+	openfpm::vector< vector_dist<3,float, aggregate<double,double>> > phases;
+
+	// first phase
+	phases.add( vector_dist<3,float, aggregate<double,double>>(0,box,bc,ghost) );
+
+	// The other 3 phases
+	phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),0) );
+	phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),0) );
+	phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),0) );
+
+	// Fill the phases with particles
+
+	auto g_it = phases.get(0).getGridIterator(sz);
+
+	while (g_it.isNext())
+	{
+		auto key = g_it.get();
+
+		// Add a particle to all the phases
+		phases.get(0).add();
+		phases.get(1).add();
+		phases.get(2).add();
+		phases.get(3).add();
+
+		phases.get(0).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(0).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(0).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
+
+		phases.get(1).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(1).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(1).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
+
+		phases.get(2).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(2).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(2).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
+
+		phases.get(3).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(3).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(3).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
+
+		++g_it;
+	}
+
+	// Sync all phases
+	for (size_t i = 0 ; i < 4 ; i++)
+	{
+		phases.get(i).map();
+		phases.get(i).ghost_get<>();
+	}
+
+	// Get the cell list of the phase 0 and 1
+	auto CL_phase0 = phases.get(0).getCellList(r_cut);
+	auto CL_phase1 = phases.get(1).getCellList(r_cut);
+
+	// This function create a Verlet-list between phases 0 and 1
+	auto NN_ver01 = createVerlet(phases.get(0),CL_phase1,r_cut);
+
+	// Check NNver0_1
+
+	bool ret = true;
+	auto it = phases.get(0).getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto p = it.get();
+		auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey());
+
+		size_t nn_count = 0;
+
+		// For each neighborhood of the particle p
+		while (Np.isNext())
+		{
+			// Neighborhood particle q
+			auto q = Np.get();
+
+			// Count the number of particles
+			nn_count++;
+
+			++Np;
+		}
+
+		ret &= nn_count == 7ul;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	// Sync all phases
+	for (size_t i = 0 ; i < 4 ; i++)
+	{
+		phases.get(i).map();
+		phases.get(i).ghost_get<>();
+	}
+
+	// NN_ver0_all
+
+	// This function create an "Empty" Multiphase Cell List
+	auto CL_all = createCellListM<2>(phases,r_cut);
+
+	// This create a Verlet-list between phase 0 and all the other phases
+	auto NNver0_all = createVerletM<2>(phases.get(0),CL_all,r_cut);
+
+	it = phases.get(0).getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto p = it.get();
+		auto Np = NNver0_all.template getNNIterator<NO_CHECK>(p.getKey());
+
+		size_t nn_cout[4] = {0,0,0,0};
+
+		// For each neighborhood of the particle p
+		while (Np.isNext())
+		{
+			// Get the particle q near to p
+			auto q = Np.getP();
+
+			// Get from which phase it come from
+			auto ph_q = Np.getV();
+
+			nn_cout[ph_q]++;
+
+			++Np;
+		}
+
+		ret &= nn_cout[0] == 7;
+		ret &= nn_cout[1] == 7;
+		ret &= nn_cout[2] == 7;
+		ret &= nn_cout[3] == 7;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+}
+
+
+BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
+{
+	if (create_vcluster().getProcessingUnits() > 24)
+		return;
+
+	size_t sz[3] = {60,60,40};
+
+	// The domain
+	Box<3,float> box({-1000.0,-1000.0,-1000.0},{2000.0,2000.0,1000.0});
+
+	// Boundary conditions
+	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+	float r_cut = 51.0;
+
+	// ghost, big enough to contain the interaction radius
+	Ghost<3,float> ghost(r_cut);
+
+	openfpm::vector< vector_dist<3,float, aggregate<size_t>> > phases;
+
+	// first phase
+	phases.add( vector_dist<3,float, aggregate<size_t>>(0,box,bc,ghost) );
+
+	// The other 3 phases
+	phases.add( vector_dist<3,float, aggregate<size_t>>(phases.get(0).getDecomposition(),0) );
+	phases.add( vector_dist<3,float, aggregate<size_t>>(phases.get(0).getDecomposition(),0) );
+	phases.add( vector_dist<3,float, aggregate<size_t>>(phases.get(0).getDecomposition(),0) );
+
+	// Fill the phases with particles
+
+	auto g_it = phases.get(0).getGridIterator(sz);
+
+	while (g_it.isNext())
+	{
+		auto key = g_it.get();
+
+		// Add a particle to all the phases
+		phases.get(0).add();
+		phases.get(1).add();
+		phases.get(2).add();
+		phases.get(3).add();
+
+		phases.get(0).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(0).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(0).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
+
+		phases.get(1).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(1).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(1).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
+
+		phases.get(2).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(2).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(2).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
+
+		phases.get(3).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(3).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(3).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
+
+		++g_it;
+	}
+
+	// Sync all phases
+	for (size_t i = 0 ; i < 4 ; i++)
+	{
+		phases.get(i).map();
+		phases.get(i).ghost_get<>();
+	}
+
+	// Get the cell list of the phase 0 and 1
+	auto CL_phase0 = phases.get(0).getCellListSym(r_cut);
+	auto CL_phase1 = phases.get(1).getCellListSym(r_cut);
+
+	// This function create a Verlet-list between phases 0 and 1
+	auto NN_ver01 = createVerletSym(phases.get(0),CL_phase1,r_cut);
+
+	// Check NNver0_1
+
+	bool ret = true;
+	auto it = phases.get(0).getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto p = it.get();
+		auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey());
+
+		// For each neighborhood of the particle p
+		while (Np.isNext())
+		{
+			// Neighborhood particle q
+			auto q = Np.get();
+
+			phases.get(0).getProp<0>(p)++;
+			phases.get(1).getProp<0>(q)++;
+
+			++Np;
+		}
+
+		++it;
+	}
+
+	phases.get(0).ghost_put<add_,0>();
+	phases.get(1).ghost_put<add_,0>();
+
+	it = phases.get(0).getDomainIterator();
+	while (it.isNext())
+	{
+		auto p = it.get();
+
+		ret &= phases.get(0).getProp<0>(p) == 7;
+		ret &= phases.get(1).getProp<0>(p) == 7;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	// Sync all phases
+	for (size_t i = 0 ; i < 4 ; i++)
+	{
+		phases.get(i).map();
+		phases.get(i).ghost_get<>();
+	}
+
+	// Reset counter on all phases
+
+	for (size_t i = 0 ; i < phases.size() ; i++)
+	{
+		it = phases.get(i).getDomainAndGhostIterator();
+		while (it.isNext())
+		{
+			auto p = it.get();
+
+			phases.get(i).getProp<0>(p) = 0;
+
+			++it;
+		}
+	}
+
+	// NN_ver0_all
+
+	// This function create an "Empty" Multiphase Cell List
+	auto CL_all = createCellListSymM<2>(phases,r_cut);
+
+	typedef decltype(createVerletSymM<2>(phases.get(0),CL_all,r_cut)) verlet_type;
+
+	verlet_type NNver_all[4];
+
+	// This create a Verlet-list between phase all phases to all the other phases
+	NNver_all[0] = createVerletSymM<2>(phases.get(0),CL_all,r_cut);
+	NNver_all[1] = createVerletSymM<2>(phases.get(1),CL_all,r_cut);
+	NNver_all[2] = createVerletSymM<2>(phases.get(2),CL_all,r_cut);
+	NNver_all[3] = createVerletSymM<2>(phases.get(3),CL_all,r_cut);
+
+	// all phases to all phases
+
+	for (size_t i = 0 ; i < phases.size() ; i++)
+	{
+		it = phases.get(i).getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto p = it.get();
+			auto Np = NNver_all[i].template getNNIterator<NO_CHECK>(p.getKey());
+
+			// For each neighborhood of the particle p
+			while (Np.isNext())
+			{
+				// Get the particle q near to p
+				auto q = Np.getP();
+
+				// Get from which phase it come from
+				auto ph_q = Np.getV();
+
+				phases.get(i).getProp<0>(p)++;
+				phases.get(ph_q).getProp<0>(q)++;
+
+				++Np;
+			}
+
+			++it;
+		}
+	}
+
+	phases.get(0).ghost_put<add_,0>();
+	phases.get(1).ghost_put<add_,0>();
+	phases.get(2).ghost_put<add_,0>();
+	phases.get(3).ghost_put<add_,0>();
+
+	it = phases.get(0).getDomainIterator();
+	while (it.isNext())
+	{
+		auto p = it.get();
+
+		ret &= phases.get(0).getProp<0>(p) == 32;
+		ret &= phases.get(1).getProp<0>(p) == 32;
+		ret &= phases.get(2).getProp<0>(p) == 32;
+		ret &= phases.get(3).getProp<0>(p) == 32;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+#endif /* SRC_VECTOR_VECTOR_DIST_MP_UNIT_TESTS_HPP_ */
diff --git a/src/Vector/vector_dist_multiphase_functions.hpp b/src/Vector/vector_dist_multiphase_functions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4a3c6e32bd51a19a6d4fddb68b1f75d398401963
--- /dev/null
+++ b/src/Vector/vector_dist_multiphase_functions.hpp
@@ -0,0 +1,126 @@
+/*
+ * vector_dist_multiphase_functions.hpp
+ *
+ *  Created on: Oct 14, 2016
+ *      Author: i-bird
+ */
+
+#ifndef SRC_VECTOR_VECTOR_DIST_MULTIPHASE_FUNCTIONS_HPP_
+#define SRC_VECTOR_VECTOR_DIST_MULTIPHASE_FUNCTIONS_HPP_
+
+#include "NN/CellList/CellListM.hpp"
+#include "NN/VerletList/VerletListM.hpp"
+
+template<typename Vector,typename CL, typename T> VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> createVerlet(Vector & v, CL & cl, T r_cut)
+{
+	VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> ver;
+
+	ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local());
+
+	return ver;
+}
+
+template<unsigned int sh_byte, typename Vector,typename CL, typename T> VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> createVerletM(Vector & v, CL & cl, T r_cut)
+{
+	VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> ver;
+
+	ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local());
+
+	return ver;
+}
+
+template<unsigned int nbit, typename Vector, typename T> CellListM<Vector::dims,typename Vector::stype,nbit> createCellListM(openfpm::vector<Vector> & phases, T r_cut)
+{
+	size_t div[3];
+	Box<Vector::dims,typename Vector::stype> box_cl;
+
+	CellListM<Vector::dims,typename Vector::stype,nbit> NN;
+	if (phases.size() == 0)
+		return NN;
+
+	box_cl = phases.get(0).getDecomposition().getProcessorBounds();
+	phases.get(0).getCellListParams(r_cut,div,box_cl);
+
+	NN.Initialize(box_cl,div);
+
+	// for all the phases i
+	for (size_t i = 0; i < phases.size() ; i++)
+	{
+		// iterate across all the particle of the phase i
+		auto it = phases.get(i).getDomainAndGhostIterator();
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			// Add the particle of the phase i to the cell list
+			NN.add(phases.get(i).getPos(key), key.getKey(), i);
+
+			++it;
+		}
+	}
+
+	return NN;
+}
+
+
+/////// Symmetric version
+
+template<typename Vector,typename CL, typename T> VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> createVerletSym(Vector & v, CL & cl, T r_cut)
+{
+	VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> ver;
+
+	ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local());
+
+	return ver;
+}
+
+template<unsigned int sh_byte, typename Vector,typename CL, typename T> VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> createVerletSymM(Vector & v, CL & cl, T r_cut)
+{
+	VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> ver;
+
+	ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local(),VL_SYMMETRIC);
+
+	return ver;
+}
+
+template<unsigned int nbit, typename Vector, typename T> CellListM<Vector::dims,typename Vector::stype,nbit> createCellListSymM(openfpm::vector<Vector> & phases, T r_cut)
+{
+	Box<Vector::dims,typename Vector::stype> box_cl;
+
+	CellListM<Vector::dims,typename Vector::stype,nbit> NN;
+	if (phases.size() == 0)
+		return NN;
+
+	// Calculate the Cell list division for this CellList
+	CellDecomposer_sm<Vector::dims,typename Vector::stype,shift<Vector::dims,typename Vector::stype>> cd_sm;
+
+	size_t pad = 0;
+	cl_param_calculateSym(phases.get(0).getDecomposition().getDomain(),cd_sm,phases.get(0).getDecomposition().getGhost(),r_cut,pad);
+
+	// Processor bounding box
+	Box<Vector::dims, typename Vector::stype> pbox = phases.get(0).getDecomposition().getProcessorBounds();
+
+	// Ghost padding extension
+	Ghost<Vector::dims,size_t> g_ext(0);
+	NN.Initialize(cd_sm,pbox,pad);
+
+	// for all the phases i
+	for (size_t i = 0; i < phases.size() ; i++)
+	{
+		// iterate across all the particle of the phase i
+		auto it = phases.get(i).getDomainAndGhostIterator();
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			// Add the particle of the phase i to the cell list
+			NN.add(phases.get(i).getPos(key), key.getKey(), i);
+
+			++it;
+		}
+	}
+
+	return NN;
+}
+
+#endif /* SRC_VECTOR_VECTOR_DIST_MULTIPHASE_FUNCTIONS_HPP_ */
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index 8c3c800e3d974914dfa011fed44c492824116f37..866868b36a6b3d2f285a76eccc23ff812164ca92 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -1523,6 +1523,24 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
 	vd.map();
 	vd.ghost_get<0,1,2>();
 
+	// shift the particle position by 1.0
+
+	it = vd.getGhostIterator();
+	while (it.isNext())
+	{
+		// Particle p
+		auto p = it.get();
+
+		// we shift down he particles
+		vd.getPos(p)[0] = 10.0;
+
+		// we shift
+		vd.getPos(p)[1] = 17.0;
+
+		// next particle
+		++it;
+	}
+
 	for (size_t i = 0 ; i < 10 ; i++)
 	{
 		auto it = vd.getDomainIterator();
@@ -1543,7 +1561,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
 			++it;
 		}
 
-		vd.ghost_get<0>(SKIP_LABELLING);
+		vd.ghost_get<0>(SKIP_LABELLING | NO_POSITION);
 
 		auto it2 = vd.getGhostIterator();
 		bool ret = true;
@@ -1560,6 +1578,22 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
 			++it2;
 		}
 
+		it2 = vd.getGhostIterator();
+		while (it2.isNext())
+		{
+			// Particle p
+			auto p = it.get();
+
+			// we shift down he particles
+			ret &= vd.getPos(p)[0] == 10.0;
+
+			// we shift
+			ret &= vd.getPos(p)[1] == 17.0;
+
+			// next particle
+			++it2;
+		}
+
 		BOOST_REQUIRE_EQUAL(ret,true);
 	}
 }