diff --git a/openfpm_numerics b/openfpm_numerics
index 9a13102746a2160146c84318e54832c4d1c8e30e..4cbdaaa7e37d236b598c130187170559d07d94c3 160000
--- a/openfpm_numerics
+++ b/openfpm_numerics
@@ -1 +1 @@
-Subproject commit 9a13102746a2160146c84318e54832c4d1c8e30e
+Subproject commit 4cbdaaa7e37d236b598c130187170559d07d94c3
diff --git a/src/Grid/staggered_dist_grid.hpp b/src/Grid/staggered_dist_grid.hpp
index 3b90cad3a65e2ba39e93c57eefe7c0f15e9c75bf..7c0a7e875d752ce0f7b833512245bd29d032f2bd 100644
--- a/src/Grid/staggered_dist_grid.hpp
+++ b/src/Grid/staggered_dist_grid.hpp
@@ -117,7 +117,7 @@ public:
 		if (extends< typename boost::mpl::at<typename T::type,boost::mpl::int_<p> >::type >::mul() != cmb.size())
 			std::cerr << __FILE__ << ":" << __LINE__ << " error properties has " << extends< typename boost::mpl::at<typename T::type,boost::mpl::int_<p> >::type >::mul() << " components, but " << cmb.size() << "has been defined \n";
 #endif
-		c_prp.get(p) = cmb;
+		c_prp[p] = cmb;
 	}
 
 	/*! \brief It set all the properties defined to be staggered on the default location
diff --git a/src/Vector/se_class3_vector.hpp b/src/Vector/se_class3_vector.hpp
index 11611166c1f083341255b5f6014e6fc2dde82d4e..2b6adae93b061d7c9e9e2a3ec9d0a7cecef198c7 100644
--- a/src/Vector/se_class3_vector.hpp
+++ b/src/Vector/se_class3_vector.hpp
@@ -16,6 +16,7 @@
 
 #define SE3_STATUS -2
 #define SE3_TYPE -1
+#define SE3_SIZE 2
 
 enum statuses
 {
@@ -37,6 +38,67 @@ enum ptype
 	INSIDE
 };
 
+// is initialized
+template<typename T>
+struct is_initialized
+{
+	static const int init = UNINITIALIZED;
+};
+
+// is initialized
+template<typename T>
+struct is_initialized<openfpm::vector<T>>
+{
+	static const int init = CLEAN;
+};
+
+
+///////////////////////////////
+
+/*! \brief this class is a functor for "for_each" algorithm
+ *
+ * This class is a functor for "for_each" algorithm. For each
+ * element of the boost::vector the operator() is called.
+ * Is mainly used to initialize the properties
+ *
+ * \tparam encap source
+ * \tparam encap dst
+ *
+ */
+
+template<unsigned int Np, typename vector>
+struct init_prop
+{
+	//! vector for prop initializetion
+	size_t (& prp_init)[Np];
+
+	/*! \brief constructor
+	 *
+	 *
+	 * \param src encapsulated object1
+	 * \param dst encapsulated object2
+	 *
+	 */
+	inline init_prop(size_t ( & prp_init)[Np])
+	:prp_init(prp_init)
+	{
+	};
+
+
+	/*!  \brief It call the copy function for each property
+	 *
+	 * \param t each member
+	 *
+	 */
+	template<typename T>
+	inline void operator()(T& t) const
+	{
+		typedef typename boost::mpl::at<vector,boost::mpl::int_<T::value>>::type tc;
+
+		prp_init[T::value] = is_initialized<tc>::init;
+	}
+};
+
 // Unknown type
 template<typename tcheck, bool foundamental>
 struct typeCheck
@@ -147,7 +209,7 @@ template<typename vector>
 struct propCheckNAN
 {
 	//! Data to check
-	vector & data;
+	const vector & data;
 
 	//! Element to check
 	size_t id;
@@ -157,7 +219,7 @@ struct propCheckNAN
 	 * \param
 	 *
 	 */
-	inline propCheckNAN(vector & data, size_t id)
+	inline propCheckNAN(const vector & data, size_t id)
 	:data(data),id(id)
 	{};
 
@@ -196,7 +258,7 @@ template<typename vector>
 struct propCheckINF
 {
 	//! Data to check
-	vector & data;
+	const vector & data;
 
 	//! id
 	size_t id;
@@ -207,7 +269,7 @@ struct propCheckINF
 	 * \param
 	 *
 	 */
-	inline propCheckINF(vector & data, size_t id)
+	inline propCheckINF(const vector & data, size_t id)
 	:data(data),id(id)
 	{};
 
@@ -261,7 +323,7 @@ class se_class3_vector
 		int sync[2][Np];
 
 		//! number of real properties + POSITION
-		static const size_t Np_real = Np+SE3_STATUS+1;
+		static const size_t Np_real = Np+SE3_STATUS;
 
 		//! Domain decomposition object
 		Decomposition & dec;
@@ -272,6 +334,9 @@ class se_class3_vector
 		//! temporal buffer
 		openfpm::vector<size_t> non_NP;
 
+		//! last write
+		size_t l_wrt;
+
 		bool isLocalHalo(const Point<dim,T> & p)
 		{
 			for (size_t i = 0; i < dec.getNLocalSub(); i++)
@@ -335,7 +400,7 @@ class se_class3_vector
 		{
 			non_NP.clear();
 
-			for (size_t i = 0 ; i < Np_real-1 ; i++)
+			for (size_t i = 0 ; i < Np_real ; i++)
 			{
 				bool found = false;
 				for (size_t j = 0 ; j < sizeof...(prp) ; j++)
@@ -355,7 +420,7 @@ class se_class3_vector
 
 		std::string getPrpName(size_t i) const
 		{
-			if (i == Np_real-1)
+			if (i == Np_real)
 				return std::string("POSITION");
 
 			return std::to_string(i);
@@ -382,8 +447,9 @@ class se_class3_vector
 			{
 				auto p = it.get();
 
-				for (size_t i = 0 ; i < Np_real ; i++)
-					vd.template getPropNC<Np+SE3_STATUS>(p)[i] = UNINITIALIZED;
+				init_prop<Np_real+1,typename vector::value_type::type> np_r(vd.template getPropNC<Np+SE3_STATUS>(p));
+
+				boost::mpl::for_each_ref< boost::mpl::range_c<int,0,Np_real+1> >(np_r);
 
 				vd.template getPropNC<Np+SE3_TYPE>(p) = INSIDE;
 
@@ -452,9 +518,12 @@ class se_class3_vector
 				for (size_t i = 0 ; i < sizeof...(prp) ; i++)
 				{
 					if (vd.template getPropNC<Np+SE3_STATUS>(p)[gg[i]] == DIRTY)
-						vd.template getPropNC<Np+SE3_STATUS>(p)[gg[i]] = CLEAN;
+					{vd.template getPropNC<Np+SE3_STATUS>(p)[gg[i]] = CLEAN;}
 				}
 
+				if (vd.template getPropNC<Np+SE3_STATUS>(p)[Np_real] == DIRTY)
+				{vd.template getPropNC<Np+SE3_STATUS>(p)[Np_real] = CLEAN;}
+
 				vd.template getPropNC<Np+SE3_TYPE>(p) = GHOST;
 
 				++it2;
@@ -463,7 +532,7 @@ class se_class3_vector
 			if (!(opt & KEEP_PROPERTIES))
 			{
 				for (size_t i = 0 ; i < non_NP.size() ; i++)
-					sync[GHOST][gg[i]] = NOTSYNC;
+					sync[GHOST][non_NP.get(i)] = NOTSYNC;
 
 				auto it = vd.getGhostIterator_no_se3();
 
@@ -484,7 +553,9 @@ class se_class3_vector
 				sync[GHOST][gg[i]] = SYNC;
 
 			if (!(opt & NO_POSITION))
+			{
 				sync[GHOST][Np_real] = SYNC;
+			}
 
 			if (!(opt & KEEP_PROPERTIES))
 			{
@@ -549,52 +620,12 @@ class se_class3_vector
 			{
 				auto p = it.get();
 
-				for (size_t j = 0 ; j < Np ; j++)
+				for (size_t j = 0 ; j < Np_real + 1 ; j++)
 				{
 					if (vd.template getPropNC<Np+SE3_STATUS>(p)[j] == DIRTY)
 						vd.template getPropNC<Np+SE3_STATUS>(p)[j] = CLEAN;
 				}
 
-#ifdef CHECK_FOR_POSINF
-
-				if ( std::isinf(vd.getPosNC(p)[0]) || std::isinf(vd.getPosNC(p)[1]) || std::isinf(vd.getPosNC(p)[2]) )
-				{
-					std::cerr << __FILE__ << ":" << __LINE__ << " error detected INF in position for particle p=" << p.getKey() << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
-					ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
-				}
-
-#endif
-
-#ifdef CHECKFOR_POSNAN
-
-				if ( std::isnan(vd.getPosNC(p)[0]) || std::isnan(vd.getPosNC(p)[1]) || std::isnan(vd.getPosNC(p)[2]) )
-				{
-					std::cerr << __FILE__ << ":" << __LINE__ << " error detected NAN in position for particle p=" << p.getKey() << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
-					ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
-				}
-
-#endif
-
-#ifdef CHECKFOR_PROPINF
-
-				{
-					propCheckINF<vector> checker(vd,p.getKey());
-
-					boost::mpl::for_each_ref< boost::mpl::range_c<int,0, Np_real-1 > > (checker);
-				}
-
-#endif
-
-#ifdef CHECKFOR_PROPNAN
-
-				{
-					propCheckNAN<vector> checker(vd,p.getKey());
-
-					boost::mpl::for_each_ref< boost::mpl::range_c<int,0, Np_real-1 > >(checker);
-				}
-
-#endif
-
 				++it;
 			}
 		}
@@ -607,7 +638,7 @@ class se_class3_vector
 			{
 				auto p = it.get();
 
-				for (size_t j = 0 ; j < Np ; j++)
+				for (size_t j = 0 ; j < Np_real + 1 ; j++)
 				{
 					if (vd.template getPropNC<Np+SE3_STATUS>(p)[j] == DIRTY)
 					{
@@ -621,7 +652,7 @@ class se_class3_vector
 
 		void map_post()
 		{
-			for (size_t j = 0 ; j < Np ; j++)
+			for (size_t j = 0 ; j < Np_real + 1 ; j++)
 			{
 
 				sync[GHOST][j] = NOTSYNC;
@@ -633,12 +664,9 @@ class se_class3_vector
 			{
 				auto p = it.get();
 
-				for (size_t j = 0 ; j < Np ; j++)
-				{
-					Point<vector::dims,typename vector::stype> xp = vd.getPosNC(p);
+				Point<vector::dims,typename vector::stype> xp = vd.getPosNC(p);
 
-					vd.template getPropNC<Np+SE3_TYPE>(p) = getParticleType(xp,p.getKey(),vd);
-				}
+				vd.template getPropNC<Np+SE3_TYPE>(p) = getParticleType(xp,p.getKey(),vd);
 
 				++it;
 			}
@@ -657,20 +685,20 @@ class se_class3_vector
 					str << __FILE__ << ":" << __LINE__ << " Error you are reading from the particle " << p << " of type " << type_str << " the property=" << getPrpName(prp) << ". But it result to be uninitialized" << std::endl;
 
 				// It is an error read from an uninitialized property
-				std::cerr << str.str();
+				std::cerr << str.str() << std::endl;
 				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
 			}
 
-			if (vd.template getPropNC<Np+SE3_STATUS>(p)[prp] == DIRTY)
+			if (vd.template getPropNC<Np+SE3_STATUS>(p)[prp] == DIRTY && p != l_wrt)
 			{
-				std::cerr << __FILE__ << ":" << __LINE__ << " Warning you are reading from an particle that has been changed already in the same cycle" << std::endl;
+				std::cerr << __FILE__ << ":" << __LINE__ << " Warning you are reading from a particle that has been changed already in the same cycle" << std::endl;
 				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
 			}
 
-			if (vd.template getPropNC<Np+SE3_TYPE>(p) == GHOST || vd.template getPropNC<Np+SE3_TYPE>(p) == HALO)
+			if (vd.template getPropNC<Np+SE3_TYPE>(p) == GHOST)
 			{
 				// if we read from the ghost we have to ensure that the ghost is in
-				// sync in particular that the state of the halo is CLEAR
+				// sync in particular that the state of the halo is CLEAN
 
 				if (sync[vd.template getPropNC<Np+SE3_TYPE>(p)][prp] != SYNC)
 				{
@@ -678,6 +706,47 @@ class se_class3_vector
 					ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
 				}
 			}
+
+#ifdef CHECK_FOR_POSINF
+
+			if ( std::isinf(vd.getPosNC(p)[0]) || std::isinf(vd.getPosNC(p)[1]) || std::isinf(vd.getPosNC(p)[2]) )
+			{
+				std::cerr << __FILE__ << ":" << __LINE__ << " error detected INF in position for particle p=" << p << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
+				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
+			}
+
+#endif
+
+#ifdef CHECKFOR_POSNAN
+
+			if ( std::isnan(vd.getPosNC(p)[0]) || std::isnan(vd.getPosNC(p)[1]) || std::isnan(vd.getPosNC(p)[2]) )
+			{
+				std::cerr << __FILE__ << ":" << __LINE__ << " error detected NAN in position for particle p=" << p << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
+				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
+			}
+
+#endif
+
+#ifdef CHECKFOR_PROPINF
+
+			{
+				propCheckINF<vector> checker(vd,p);
+
+				boost::mpl::for_each_ref< boost::mpl::range_c<int,0, Np_real > > (checker);
+			}
+
+#endif
+
+#ifdef CHECKFOR_PROPNAN
+
+			{
+				propCheckNAN<vector> checker(vd,p);
+
+				boost::mpl::for_each_ref< boost::mpl::range_c<int,0, Np_real > >(checker);
+			}
+
+#endif
+
 		}
 
 		template<unsigned int prp> void write(vector & vd, size_t p)
@@ -691,33 +760,29 @@ class se_class3_vector
 					vd.get_se_class3().template setGhostOutSync<prp>();
 			}
 
+			l_wrt = p;
 		}
 
 		//! Copy operator
 		se_class3_vector<Np,dim,T,Decomposition,vector> & operator=(const se_class3_vector<Np,dim,T,Decomposition,vector> & se3)
 		{
-			for (size_t i = 0 ; i < Np ; i++)
+			for (size_t i = 0 ; i < Np_real + 1 ; i++)
 			{
 				sync[0][i] = se3.sync[0][i];
 				sync[1][i] = se3.sync[1][i];
 			}
 
-			dec = se3.dec;
-			vd = se3.vd;
-
 			return *this;
 		}
 
 		template<unsigned int prp> void setHaloOutSync()
 		{
-			for (size_t i = 0 ; i < Np ; i++)
-				sync[HALO][prp] = NOTSYNC;
+			sync[HALO][prp] = NOTSYNC;
 		}
 
 		template<unsigned int prp> void setGhostOutSync()
 		{
-			for (size_t i = 0 ; i < Np ; i++)
-				sync[GHOST][prp] = NOTSYNC;
+			sync[GHOST][prp] = NOTSYNC;
 		}
 
 		void getNN()
diff --git a/src/Vector/se_class3_vector_unit_tests.hpp b/src/Vector/se_class3_vector_unit_tests.hpp
index 32d6ef55a782e7c638c35b3c9d8a69c9045b5849..de8261e5bdced1e60007620a75916758d04aba15 100644
--- a/src/Vector/se_class3_vector_unit_tests.hpp
+++ b/src/Vector/se_class3_vector_unit_tests.hpp
@@ -558,7 +558,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_class3_check_nan_inf )
 	{
 	try
 	{
-		auto it = vd.getDomainIterator();
+		vd.getPosRead(0);
 	}
 	catch (std::exception & e)
 	{
@@ -585,7 +585,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_class3_check_nan_inf )
 	{
 	try
 	{
-		auto it = vd2.getDomainIterator();
+		vd2.getLastPosRead();
 	}
 	catch (std::exception & e)
 	{
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 69ef226ed27c62587c332818a099ddbc07443484..f981d94658ecb8af331406310ea5b80c87d96624 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -743,6 +743,8 @@ public:
 #endif
 	}
 
+#ifndef ONLY_READWRITE_GETTER
+
 	/*! \brief Get the position of the last element
 	 *
 	 * \return the position of the element in space
@@ -765,6 +767,8 @@ public:
 		return v_prp.template get<id>(g_m - 1);
 	}
 
+#endif
+
 ////////////////////////////// READ AND WRITE VERSION //////////
 
 	/*! \brief Get the position of the last element
@@ -775,7 +779,7 @@ public:
 	inline auto getLastPosRead() -> decltype(v_pos.template get<0>(0))
 	{
 #ifdef SE_CLASS3
-		se3.read<prop::max_prop_real>(*this,g_m-1);
+		se3.template read<prop::max_prop_real>(*this,g_m-1);
 #endif
 
 		return v_pos.template get<0>(g_m - 1);
@@ -806,7 +810,7 @@ public:
 	inline auto getLastPosWrite() -> decltype(v_pos.template get<0>(0))
 	{
 #ifdef SE_CLASS3
-		se3.write<prop::max_prop_real>(*this,g_m-1);
+		se3.template write<prop::max_prop_real>(*this,g_m-1);
 #endif
 
 		return v_pos.template get<0>(g_m - 1);
@@ -822,7 +826,7 @@ public:
 	template<unsigned int id> inline auto getLastPropWrite() -> decltype(v_prp.template get<id>(0))
 	{
 #ifdef SE_CLASS3
-		se3.write<id>(*this,g_m-1);
+		se3.template write<id>(*this,g_m-1);
 #endif
 
 		return v_prp.template get<id>(g_m - 1);
@@ -882,10 +886,11 @@ public:
 	 * \return the Cell list
 	 *
 	 */
-	template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast, shift<dim, St> > > CellL getCellList(St r_cut)
+	template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast, shift<dim, St> > > CellL getCellList(St r_cut, bool no_se3 = false)
 	{
 #ifdef SE_CLASS3
-		se3.getNN();
+		if (no_se3 == false)
+		{se3.getNN();}
 #endif
 #ifdef SE_CLASS1
 		check_ghost_compatible_rcut(r_cut);
@@ -895,7 +900,7 @@ public:
 		Ghost<dim,St> g = getDecomposition().getGhost();
 		g.magnify(1.013);
 
-		return getCellList<CellL>(r_cut, g);
+		return getCellList<CellL>(r_cut, g,no_se3);
 	}
 
 	/*! \brief Construct an hilbert cell list starting from the stored particles
@@ -930,10 +935,11 @@ public:
 	 * \param cell_list Cell list to update
 	 *
 	 */
-	template<typename CellL> void updateCellList(CellL & cell_list)
+	template<typename CellL> void updateCellList(CellL & cell_list, bool no_se3 = false)
 	{
 #ifdef SE_CLASS3
-		se3.getNN();
+		if (no_se3 == false)
+		{se3.getNN();}
 #endif
 
 		// This function assume equal spacing in all directions
@@ -1012,10 +1018,11 @@ public:
 	 * \return the CellList
 	 *
 	 */
-	template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast, shift<dim, St> > > CellL getCellList(St r_cut, const Ghost<dim, St> & enlarge)
+	template<typename CellL = CellList_gen<dim, St, Process_keys_lin, Mem_fast, shift<dim, St> > > CellL getCellList(St r_cut, const Ghost<dim, St> & enlarge, bool no_se3 = false)
 	{
 #ifdef SE_CLASS3
-		se3.getNN();
+		if (no_se3 == false)
+		{se3.getNN();}
 #endif
 
 		CellL cell_list;
@@ -1032,7 +1039,7 @@ public:
 		cell_list.Initialize(pbox, div, g_m);
 		cell_list.set_ndec(getDecomposition().get_ndec());
 
-		updateCellList(cell_list);
+		updateCellList(cell_list,no_se3);
 
 		return cell_list;
 	}
@@ -1987,6 +1994,10 @@ public:
 	 */
 	template<typename cli> ParticleItCRS_Cells<dim,cli> getParticleIteratorCRS_Cell(cli & NN)
 	{
+#ifdef SE_CLASS3
+		se3.getIterator();
+#endif
+
 #ifdef SE_CLASS1
 		if (!(opt & BIND_DEC_TO_GHOST))
 		{
diff --git a/src/Vector/vector_dist_HDF5_chckpnt_restart_test.hpp b/src/Vector/vector_dist_HDF5_chckpnt_restart_test.hpp
index 1ed0477ec6a402a514a00508d72e55fb96446187..ca4b26ce0f33b19898da484ed18613bc43ef7b75 100644
--- a/src/Vector/vector_dist_HDF5_chckpnt_restart_test.hpp
+++ b/src/Vector/vector_dist_HDF5_chckpnt_restart_test.hpp
@@ -119,6 +119,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_hdf5_save_test )
 
 BOOST_AUTO_TEST_CASE( vector_dist_hdf5_load_test )
 {
+#ifndef SE_CLASS3
+
 	Vcluster & v_cl = create_vcluster();
 
 	Box<dim,float> box;
@@ -169,6 +171,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_hdf5_load_test )
 
 		++it2;
 	}
+
+#endif
 }
 
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/Vector/vector_dist_MP_unit_tests.hpp b/src/Vector/vector_dist_MP_unit_tests.hpp
index 74414c67ceb518c2789493ba9e9aa6dba0b1150f..807b7a259f435db2839f345097e2e343d38542f8 100644
--- a/src/Vector/vector_dist_MP_unit_tests.hpp
+++ b/src/Vector/vector_dist_MP_unit_tests.hpp
@@ -222,21 +222,21 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
 		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(0).getLastPosWrite()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(0).getLastPosWrite()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(0).getLastPosWrite()[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(1).getLastPosWrite()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(1).getLastPosWrite()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(1).getLastPosWrite()[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(2).getLastPosWrite()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(2).getLastPosWrite()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(2).getLastPosWrite()[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);
+		phases.get(3).getLastPosWrite()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
+		phases.get(3).getLastPosWrite()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
+		phases.get(3).getLastPosWrite()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
 
 		++g_it;
 	}
@@ -289,8 +289,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
 			// Neighborhood particle q
 			auto q = Np.get();
 
-			phases.get(0).getProp<0>(p)++;
-			phases.get(1).getProp<0>(q)++;
+			phases.get(0).getPropWrite<0>(p)++;
+			phases.get(1).getPropWrite<0>(q)++;
 
 			++Np;
 		}
@@ -301,13 +301,19 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
 	phases.get(0).ghost_put<add_,0>();
 	phases.get(1).ghost_put<add_,0>();
 
+#ifdef SE_CLASS3
+
+	phases.get(1).getDomainIterator();
+
+#endif
+
 	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;
+		ret &= phases.get(0).getPropRead<0>(p) == 7;
+		ret &= phases.get(1).getPropRead<0>(p) == 7;
 
 		++it;
 	}
@@ -330,7 +336,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
 		{
 			auto p = it.get();
 
-			phases.get(i).getProp<0>(p) = 0;
+			phases.get(i).getPropWrite<0>(p) = 0;
 
 			++it;
 		}
@@ -371,8 +377,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
 				// 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)++;
+				phases.get(i).getPropWrite<0>(p)++;
+				phases.get(ph_q).getPropWrite<0>(q)++;
 
 				++Np;
 			}
@@ -386,22 +392,30 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
 	phases.get(2).ghost_put<add_,0>();
 	phases.get(3).ghost_put<add_,0>();
 
+#ifdef SE_CLASS3
+
+	it = phases.get(1).getDomainIterator();
+	it = phases.get(2).getDomainIterator();
+	it = phases.get(3).getDomainIterator();
+
+#endif
+
 	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;
+		ret &= phases.get(0).getPropRead<0>(p) == 32;
+		ret &= phases.get(1).getPropRead<0>(p) == 32;
+		ret &= phases.get(2).getPropRead<0>(p) == 32;
+		ret &= phases.get(3).getPropRead<0>(p) == 32;
 
 		if (ret == false)
 		{
-			std::cout << phases.get(0).getProp<0>(p) << std::endl;
-			std::cout << phases.get(1).getProp<0>(p) << std::endl;
-			std::cout << phases.get(2).getProp<0>(p) << std::endl;
-			std::cout << phases.get(3).getProp<0>(p) << std::endl;
+			std::cout << phases.get(0).getPropRead<0>(p) << std::endl;
+			std::cout << phases.get(1).getPropRead<0>(p) << std::endl;
+			std::cout << phases.get(2).getPropRead<0>(p) << std::endl;
+			std::cout << phases.get(3).getPropRead<0>(p) << std::endl;
 		}
 
 		++it;
diff --git a/src/Vector/vector_dist_cell_list_tests.hpp b/src/Vector/vector_dist_cell_list_tests.hpp
index 7603d87d91bd749552b3d4963c5c1e1c50907314..0962a3cd5f2439fb252c01aceefe21eaeeb195a8 100644
--- a/src/Vector/vector_dist_cell_list_tests.hpp
+++ b/src/Vector/vector_dist_cell_list_tests.hpp
@@ -63,9 +63,11 @@ BOOST_AUTO_TEST_CASE( vector_dist_reorder_2d_test )
 
 		vd.map();
 
+		// in case of SE_CLASS3 get a cell-list without ghost get is an error
+
 		// Create first cell list
 
-		auto NN1 = vd.getCellList(0.01);
+		auto NN1 = vd.getCellList(0.01,true);
 
 		//An order of a curve
 		int32_t m = 6;
@@ -74,7 +76,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_reorder_2d_test )
 		vd.reorder(m);
 
 		// Create second cell list
-		auto NN2 = vd.getCellList(0.01);
+		auto NN2 = vd.getCellList(0.01,true);
 
 		//Check equality of cell sizes
 		for (size_t i = 0 ; i < NN1.getGrid().size() ; i++)
@@ -395,15 +397,15 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 	{
 		auto key = it.get();
 
-		vd.getPos(key)[0] = ud(eg);
-		vd.getPos(key)[1] = ud(eg);
-		vd.getPos(key)[2] = ud(eg);
+		vd.getPosWrite(key)[0] = ud(eg);
+		vd.getPosWrite(key)[1] = ud(eg);
+		vd.getPosWrite(key)[2] = ud(eg);
 
 		// Fill some properties randomly
 
-		vd.getProp<0>(key) = 0;
-		vd.getProp<1>(key) = 0;
-		vd.getProp<2>(key) = key.getKey() + start;
+		vd.getPropWrite<0>(key) = 0;
+		vd.getPropWrite<1>(key) = 0;
+		vd.getPropWrite<2>(key) = key.getKey() + start;
 
 		++it;
 	}
@@ -420,7 +422,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 	{
 		auto p = p_it.get();
 
-		Point<3,float> xp = vd.getPos(p);
+		Point<3,float> xp = vd.getPosRead(p);
 
 		auto Np = NN.getNNIterator(NN.getCell(xp));
 
@@ -436,7 +438,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 
 			// repulsive
 
-			Point<3,float> xq = vd.getPos(q);
+			Point<3,float> xq = vd.getPosRead(q);
 			Point<3,float> f = (xp - xq);
 
 			float distance = f.norm();
@@ -445,10 +447,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 
 			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);
+				vd.getPropWrite<0>(p)++;
+				vd.getPropWrite<3>(p).add();
+				vd.getPropWrite<3>(p).last().xq = xq;
+				vd.getPropWrite<3>(p).last().id = vd.getPropWrite<2>(q);
 			}
 
 			++Np;
@@ -467,7 +469,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 	{
 		auto p = p_it2.get();
 
-		Point<3,float> xp = vd.getPos(p);
+		Point<3,float> xp = vd.getPosRead(p);
 
 		auto Np = NN2.getNNIteratorSym<NO_CHECK>(NN2.getCell(xp),p.getKey(),vd.getPosVector());
 
@@ -483,7 +485,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 
 			// repulsive
 
-			Point<3,float> xq = vd.getPos(q);
+			Point<3,float> xq = vd.getPosRead(q);
 			Point<3,float> f = (xp - xq);
 
 			float distance = f.norm();
@@ -492,16 +494,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 
 			if (distance < r_cut )
 			{
-				vd.getProp<1>(p)++;
-				vd.getProp<1>(q)++;
+				vd.getPropWrite<1>(p)++;
+				vd.getPropWrite<1>(q)++;
 
-				vd.getProp<4>(p).add();
-				vd.getProp<4>(q).add();
+				vd.getPropWrite<4>(p).add();
+				vd.getPropWrite<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);
+				vd.getPropWrite<4>(p).last().xq = xq;
+				vd.getPropWrite<4>(q).last().xq = xp;
+				vd.getPropWrite<4>(p).last().id = vd.getProp<2>(q);
+				vd.getPropWrite<4>(q).last().id = vd.getProp<2>(p);
 			}
 
 			++Np;
@@ -520,24 +522,24 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 	{
 		auto p = p_it3.get();
 
-		ret &= vd.getProp<1>(p) == vd.getProp<0>(p);
+		ret &= vd.getPropRead<1>(p) == vd.getPropRead<0>(p);
 
-		vd.getProp<3>(p).sort();
-		vd.getProp<4>(p).sort();
+		vd.getPropRead<3>(p).sort();
+		vd.getPropRead<4>(p).sort();
 
-		ret &= vd.getProp<3>(p).size() == vd.getProp<4>(p).size();
+		ret &= vd.getPropRead<3>(p).size() == vd.getPropRead<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;
+		for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
+			ret &= vd.getPropRead<3>(p).get(i).id == vd.getPropRead<4>(p).get(i).id;
 
 		if (ret == false)
 		{
-			std::cout << vd.getProp<3>(p).size() << "   " << vd.getProp<4>(p).size() << std::endl;
+			std::cout << vd.getPropRead<3>(p).size() << "   " << vd.getPropRead<4>(p).size() << std::endl;
 
-			for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++)
-				std::cout << vd.getProp<3>(p).get(i).id << "    " << vd.getProp<4>(p).get(i).id << std::endl;
+			for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
+				std::cout << vd.getPropRead<3>(p).get(i).id << "    " << vd.getPropRead<4>(p).get(i).id << std::endl;
 
-			std::cout << vd.getProp<1>(p) << "  A  " << vd.getProp<0>(p) << std::endl;
+			std::cout << vd.getPropRead<1>(p) << "  A  " << vd.getPropRead<0>(p) << std::endl;
 
 			break;
 		}
@@ -609,23 +611,23 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
 	{
 		auto key = it.get();
 
-		vd.getPos(key)[0] = ud(eg);
-		vd.getPos(key)[1] = ud(eg);
-		vd.getPos(key)[2] = ud(eg);
+		vd.getPosWrite(key)[0] = ud(eg);
+		vd.getPosWrite(key)[1] = ud(eg);
+		vd.getPosWrite(key)[2] = ud(eg);
 
-		vd2.getPos(key)[0] = vd.getPos(key)[0];
-		vd2.getPos(key)[1] = vd.getPos(key)[1];
-		vd2.getPos(key)[2] = vd.getPos(key)[2];
+		vd2.getPosWrite(key)[0] = vd.getPosRead(key)[0];
+		vd2.getPosWrite(key)[1] = vd.getPosRead(key)[1];
+		vd2.getPosWrite(key)[2] = vd.getPosRead(key)[2];
 
 		// Fill some properties randomly
 
-		vd.getProp<0>(key) = 0;
-		vd.getProp<1>(key) = 0;
-		vd.getProp<2>(key) = key.getKey() + start;
+		vd.getPropWrite<0>(key) = 0;
+		vd.getPropWrite<1>(key) = 0;
+		vd.getPropWrite<2>(key) = key.getKey() + start;
 
-		vd2.getProp<0>(key) = 0;
-		vd2.getProp<1>(key) = 0;
-		vd2.getProp<2>(key) = key.getKey() + start;
+		vd2.getPropWrite<0>(key) = 0;
+		vd2.getPropWrite<1>(key) = 0;
+		vd2.getPropWrite<2>(key) = key.getKey() + start;
 
 		++it;
 	}
@@ -647,7 +649,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
 	{
 		auto p = p_it.get();
 
-		Point<3,float> xp = vd.getPos(p);
+		Point<3,float> xp = vd.getPosRead(p);
 
 		auto Np = NN.getNNIterator(NN.getCell(xp));
 
@@ -663,7 +665,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
 
 			// repulsive
 
-			Point<3,float> xq = vd.getPos(q);
+			Point<3,float> xq = vd.getPosRead(q);
 			Point<3,float> f = (xp - xq);
 
 			float distance = f.norm();
@@ -672,10 +674,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
 
 			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);
+				vd.getPropWrite<0>(p)++;
+				vd.getPropWrite<3>(p).add();
+				vd.getPropWrite<3>(p).last().xq = xq;
+				vd.getPropWrite<3>(p).last().id = vd.getPropRead<2>(q);
 			}
 
 			++Np;
@@ -697,7 +699,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
 	{
 		auto p = p_it2.get();
 
-		Point<3,float> xp = vd2.getPos(p);
+		Point<3,float> xp = vd2.getPosRead(p);
 
 		auto Np = p_it2.getNNIteratorCSR(vd2.getPosVector());
 
@@ -713,7 +715,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
 
 			// repulsive
 
-			Point<3,float> xq = vd2.getPos(q);
+			Point<3,float> xq = vd2.getPosRead(q);
 			Point<3,float> f = (xp - xq);
 
 			float distance = f.norm();
@@ -722,16 +724,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
 
 			if (distance < r_cut )
 			{
-				vd2.getProp<1>(p)++;
-				vd2.getProp<1>(q)++;
+				vd2.getPropWrite<1>(p)++;
+				vd2.getPropWrite<1>(q)++;
 
-				vd2.getProp<4>(p).add();
-				vd2.getProp<4>(q).add();
+				vd2.getPropWrite<4>(p).add();
+				vd2.getPropWrite<4>(q).add();
 
-				vd2.getProp<4>(p).last().xq = xq;
-				vd2.getProp<4>(q).last().xq = xp;
-				vd2.getProp<4>(p).last().id = vd2.getProp<2>(q);
-				vd2.getProp<4>(q).last().id = vd2.getProp<2>(p);
+				vd2.getPropWrite<4>(p).last().xq = xq;
+				vd2.getPropWrite<4>(q).last().xq = xp;
+				vd2.getPropWrite<4>(p).last().id = vd2.getPropRead<2>(q);
+				vd2.getPropWrite<4>(q).last().id = vd2.getPropRead<2>(p);
 			}
 
 			++Np;
@@ -743,6 +745,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
 	vd2.ghost_put<add_,1>(NO_CHANGE_ELEMENTS);
 	vd2.ghost_put<merge_,4>();
 
+#ifdef SE_CLASS3
+	vd2.getDomainIterator();
+#endif
+
 	auto p_it3 = vd.getDomainIterator();
 
 	bool ret = true;
@@ -750,16 +756,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
 	{
 		auto p = p_it3.get();
 
-		ret &= vd2.getProp<1>(p) == vd.getProp<0>(p);
+		ret &= vd2.getPropRead<1>(p) == vd.getPropRead<0>(p);
 
 
-		vd.getProp<3>(p).sort();
-		vd2.getProp<4>(p).sort();
+		vd.getPropWrite<3>(p).sort();
+		vd2.getPropWrite<4>(p).sort();
 
-		ret &= vd.getProp<3>(p).size() == vd2.getProp<4>(p).size();
+		ret &= vd.getPropRead<3>(p).size() == vd2.getPropRead<4>(p).size();
 
-		for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++)
-			ret &= vd.getProp<3>(p).get(i).id == vd2.getProp<4>(p).get(i).id;
+		for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
+			ret &= vd.getPropRead<3>(p).get(i).id == vd2.getPropRead<4>(p).get(i).id;
 
 		if (ret == false)
 			break;
@@ -827,15 +833,15 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
 	{
 		auto key = it.get();
 
-		vd.getPos(key)[0] = ud(eg);
-		vd.getPos(key)[1] = ud(eg);
-		vd.getPos(key)[2] = ud(eg);
+		vd.getPosWrite(key)[0] = ud(eg);
+		vd.getPosWrite(key)[1] = ud(eg);
+		vd.getPosWrite(key)[2] = ud(eg);
 
 		// Fill some properties randomly
 
-		vd.getProp<0>(key) = 0;
-		vd.getProp<1>(key) = 0;
-		vd.getProp<2>(key) = key.getKey() + start;
+		vd.getPropWrite<0>(key) = 0;
+		vd.getPropWrite<1>(key) = 0;
+		vd.getPropWrite<2>(key) = key.getKey() + start;
 
 		++it;
 	}
@@ -852,7 +858,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
 	{
 		auto p = p_it.get();
 
-		Point<3,float> xp = vd.getPos(p);
+		Point<3,float> xp = vd.getPosRead(p);
 
 		auto Np = NN.getNNIterator(p.getKey());
 
@@ -868,7 +874,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
 
 			// repulsive
 
-			Point<3,float> xq = vd.getPos(q);
+			Point<3,float> xq = vd.getPosRead(q);
 			Point<3,float> f = (xp - xq);
 
 			float distance = f.norm();
@@ -877,10 +883,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
 
 			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);
+				vd.getPropWrite<0>(p)++;
+				vd.getPropWrite<3>(p).add();
+				vd.getPropWrite<3>(p).last().xq = xq;
+				vd.getPropWrite<3>(p).last().id = vd.getPropRead<2>(q);
 			}
 
 			++Np;
@@ -899,7 +905,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
 	{
 		auto p = p_it2.get();
 
-		Point<3,float> xp = vd.getPos(p);
+		Point<3,float> xp = vd.getPosRead(p);
 
 		auto Np = NN2.getNNIterator<NO_CHECK>(p.getKey());
 
@@ -915,7 +921,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
 
 			// repulsive
 
-			Point<3,float> xq = vd.getPos(q);
+			Point<3,float> xq = vd.getPosRead(q);
 			Point<3,float> f = (xp - xq);
 
 			float distance = f.norm();
@@ -924,16 +930,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
 
 			if (distance < r_cut )
 			{
-				vd.getProp<1>(p)++;
-				vd.getProp<1>(q)++;
+				vd.getPropWrite<1>(p)++;
+				vd.getPropWrite<1>(q)++;
 
-				vd.getProp<4>(p).add();
-				vd.getProp<4>(q).add();
+				vd.getPropWrite<4>(p).add();
+				vd.getPropWrite<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);
+				vd.getPropWrite<4>(p).last().xq = xq;
+				vd.getPropWrite<4>(q).last().xq = xp;
+				vd.getPropWrite<4>(p).last().id = vd.getPropRead<2>(q);
+				vd.getPropWrite<4>(q).last().id = vd.getPropRead<2>(p);
 			}
 
 			++Np;
@@ -952,15 +958,15 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
 	{
 		auto p = p_it3.get();
 
-		ret &= vd.getProp<1>(p) == vd.getProp<0>(p);
+		ret &= vd.getPropRead<1>(p) == vd.getPropRead<0>(p);
 
-		vd.getProp<3>(p).sort();
-		vd.getProp<4>(p).sort();
+		vd.getPropWrite<3>(p).sort();
+		vd.getPropWrite<4>(p).sort();
 
-		ret &= vd.getProp<3>(p).size() == vd.getProp<4>(p).size();
+		ret &= vd.getPropRead<3>(p).size() == vd.getPropRead<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;
+		for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
+			ret &= vd.getPropRead<3>(p).get(i).id == vd.getPropRead<4>(p).get(i).id;
 
 		if (ret == false)
 			break;
@@ -1035,23 +1041,23 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 		{
 			auto key = it.get();
 
-			vd.getPos(key)[0] = ud(eg);
-			vd.getPos(key)[1] = ud(eg);
-			vd.getPos(key)[2] = ud(eg);
+			vd.getPosWrite(key)[0] = ud(eg);
+			vd.getPosWrite(key)[1] = ud(eg);
+			vd.getPosWrite(key)[2] = ud(eg);
 
-			vd2.getPos(key)[0] = vd.getPos(key)[0];
-			vd2.getPos(key)[1] = vd.getPos(key)[1];
-			vd2.getPos(key)[2] = vd.getPos(key)[2];
+			vd2.getPosWrite(key)[0] = vd.getPosRead(key)[0];
+			vd2.getPosWrite(key)[1] = vd.getPosRead(key)[1];
+			vd2.getPosWrite(key)[2] = vd.getPosRead(key)[2];
 
 			// Fill some properties randomly
 
-			vd.getProp<0>(key) = 0;
-			vd.getProp<1>(key) = 0;
-			vd.getProp<2>(key) = key.getKey() + start;
+			vd.getPropWrite<0>(key) = 0;
+			vd.getPropWrite<1>(key) = 0;
+			vd.getPropWrite<2>(key) = key.getKey() + start;
 
-			vd2.getProp<0>(key) = 0;
-			vd2.getProp<1>(key) = 0;
-			vd2.getProp<2>(key) = key.getKey() + start;
+			vd2.getPropWrite<0>(key) = 0;
+			vd2.getPropWrite<1>(key) = 0;
+			vd2.getPropWrite<2>(key) = key.getKey() + start;
 
 			++it;
 		}
@@ -1070,7 +1076,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 		{
 			auto p = p_it.get();
 
-			Point<3,float> xp = vd.getPos(p);
+			Point<3,float> xp = vd.getPosRead(p);
 
 			auto Np = NN.getNNIterator(p.getKey());
 
@@ -1086,7 +1092,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 
 				// repulsive
 
-				Point<3,float> xq = vd.getPos(q);
+				Point<3,float> xq = vd.getPosRead(q);
 				Point<3,float> f = (xp - xq);
 
 				float distance = f.norm();
@@ -1095,10 +1101,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 
 				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);
+					vd.getPropWrite<0>(p)++;
+					vd.getPropWrite<3>(p).add();
+					vd.getPropWrite<3>(p).last().xq = xq;
+					vd.getPropWrite<3>(p).last().id = vd.getPropRead<2>(q);
 				}
 
 				++Np;
@@ -1117,7 +1123,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 		{
 			auto p = p_it2.get();
 
-			Point<3,float> xp = vd2.getPos(p);
+			Point<3,float> xp = vd2.getPosRead(p);
 
 			auto Np = NN2.getNNIterator<NO_CHECK>(p.getKey());
 
@@ -1133,7 +1139,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 
 				// repulsive
 
-				Point<3,float> xq = vd2.getPos(q);
+				Point<3,float> xq = vd2.getPosRead(q);
 				Point<3,float> f = (xp - xq);
 
 				float distance = f.norm();
@@ -1142,16 +1148,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 
 				if (distance < r_cut )
 				{
-					vd2.getProp<1>(p)++;
-					vd2.getProp<1>(q)++;
+					vd2.getPropWrite<1>(p)++;
+					vd2.getPropWrite<1>(q)++;
 
-					vd2.getProp<4>(p).add();
-					vd2.getProp<4>(q).add();
+					vd2.getPropWrite<4>(p).add();
+					vd2.getPropWrite<4>(q).add();
 
-					vd2.getProp<4>(p).last().xq = xq;
-					vd2.getProp<4>(q).last().xq = xp;
-					vd2.getProp<4>(p).last().id = vd2.getProp<2>(q);
-					vd2.getProp<4>(q).last().id = vd2.getProp<2>(p);
+					vd2.getPropWrite<4>(p).last().xq = xq;
+					vd2.getPropWrite<4>(q).last().xq = xp;
+					vd2.getPropWrite<4>(p).last().id = vd2.getPropRead<2>(q);
+					vd2.getPropWrite<4>(q).last().id = vd2.getPropRead<2>(p);
 				}
 
 				++Np;
@@ -1164,6 +1170,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 		vd2.ghost_put<add_,1>();
 		vd2.ghost_put<merge_,4>();
 
+#ifdef SE_CLASS3
+		vd2.getDomainIterator();
+#endif
+
 		auto p_it3 = vd.getDomainIterator();
 
 		bool ret = true;
@@ -1171,16 +1181,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 		{
 			auto p = p_it3.get();
 
-			ret &= vd2.getProp<1>(p) == vd.getProp<0>(p);
+			ret &= vd2.getPropRead<1>(p) == vd.getPropRead<0>(p);
 
 
-			vd.getProp<3>(p).sort();
-			vd2.getProp<4>(p).sort();
+			vd.getPropWrite<3>(p).sort();
+			vd2.getPropWrite<4>(p).sort();
 
-			ret &= vd.getProp<3>(p).size() == vd2.getProp<4>(p).size();
+			ret &= vd.getPropRead<3>(p).size() == vd2.getPropRead<4>(p).size();
 
-			for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++)
-				ret &= vd.getProp<3>(p).get(i).id == vd2.getProp<4>(p).get(i).id;
+			for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
+				ret &= vd.getPropRead<3>(p).get(i).id == vd2.getPropRead<4>(p).get(i).id;
 
 			if (ret == false)
 				break;
@@ -1205,23 +1215,23 @@ template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop >
 	{
 		auto key = it.get();
 
-		vd.getPos(key)[0] = ud(eg);
-		vd.getPos(key)[1] = ud(eg);
-		vd.getPos(key)[2] = ud(eg);
+		vd.getPosWrite(key)[0] = ud(eg);
+		vd.getPosWrite(key)[1] = ud(eg);
+		vd.getPosWrite(key)[2] = ud(eg);
 
-		vd2.getPos(key)[0] = vd.getPos(key)[0];
-		vd2.getPos(key)[1] = vd.getPos(key)[1];
-		vd2.getPos(key)[2] = vd.getPos(key)[2];
+		vd2.getPosWrite(key)[0] = vd.getPosRead(key)[0];
+		vd2.getPosWrite(key)[1] = vd.getPosRead(key)[1];
+		vd2.getPosWrite(key)[2] = vd.getPosRead(key)[2];
 
 		// Fill some properties randomly
 
-		vd.template getProp<0>(key) = 0;
-		vd.template getProp<1>(key) = 0;
-		vd.template getProp<2>(key) = key.getKey() + start;
+		vd.template getPropWrite<0>(key) = 0;
+		vd.template getPropWrite<1>(key) = 0;
+		vd.template getPropWrite<2>(key) = key.getKey() + start;
 
-		vd2.template getProp<0>(key) = 0;
-		vd2.template getProp<1>(key) = 0;
-		vd2.template getProp<2>(key) = key.getKey() + start;
+		vd2.template getPropWrite<0>(key) = 0;
+		vd2.template getPropWrite<1>(key) = 0;
+		vd2.template getPropWrite<2>(key) = key.getKey() + start;
 
 		++it;
 	}
@@ -1229,8 +1239,6 @@ template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop >
 	vd.map();
 	vd2.map();
 
-	Vcluster & v_cl = create_vcluster();
-
 	// sync the ghost
 	vd.template ghost_get<0,2>();
 	vd2.template ghost_get<0,2>();
@@ -1242,14 +1250,7 @@ template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop >
 	{
 		auto p = p_it.get();
 
-		Point<3,float> xp = vd.getPos(p);
-
-		if (v_cl.getProcessUnitID() == 2 && p.getKey() == 137)
-		{
-			int debug = 0;
-			debug++;
-		}
-
+		Point<3,float> xp = vd.getPosRead(p);
 
 		auto Np = NN.getNNIterator(p.getKey());
 
@@ -1265,7 +1266,7 @@ template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop >
 
 			// repulsive
 
-			Point<3,float> xq = vd.getPos(q);
+			Point<3,float> xq = vd.getPosRead(q);
 			Point<3,float> f = (xp - xq);
 
 			float distance = f.norm();
@@ -1274,10 +1275,10 @@ template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop >
 
 			if (distance < r_cut )
 			{
-				vd.template getProp<0>(p)++;
-				vd.template getProp<3>(p).add();
-				vd.template getProp<3>(p).last().xq = xq;
-				vd.template getProp<3>(p).last().id = vd.template getProp<2>(q);
+				vd.template getPropWrite<0>(p)++;
+				vd.template getPropWrite<3>(p).add();
+				vd.template getPropWrite<3>(p).last().xq = xq;
+				vd.template getPropWrite<3>(p).last().id = vd.template getPropRead<2>(q);
 			}
 
 			++Np;
@@ -1297,7 +1298,7 @@ template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop >
 	{
 		auto p = p_it2.get();
 
-		Point<3,float> xp = vd2.getPos(p);
+		Point<3,float> xp = vd2.getPosRead(p);
 
 		auto Np = NN2.template getNNIterator<NO_CHECK>(p);
 
@@ -1313,23 +1314,23 @@ template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop >
 
 			// repulsive
 
-			Point<3,float> xq = vd2.getPos(q);
+			Point<3,float> xq = vd2.getPosRead(q);
 			Point<3,float> f = (xp - xq);
 
 			float distance = f.norm();
 
 			if (distance < r_cut )
 			{
-				vd2.template getProp<1>(p)++;
-				vd2.template getProp<1>(q)++;
+				vd2.template getPropWrite<1>(p)++;
+				vd2.template getPropWrite<1>(q)++;
 
-				vd2.template getProp<4>(p).add();
-				vd2.template getProp<4>(q).add();
+				vd2.template getPropWrite<4>(p).add();
+				vd2.template getPropWrite<4>(q).add();
 
-				vd2.template getProp<4>(p).last().xq = xq;
-				vd2.template getProp<4>(q).last().xq = xp;
-				vd2.template getProp<4>(p).last().id = vd2.template getProp<2>(q);
-				vd2.template getProp<4>(q).last().id = vd2.template getProp<2>(p);
+				vd2.template getPropWrite<4>(p).last().xq = xq;
+				vd2.template getPropWrite<4>(q).last().xq = xp;
+				vd2.template getPropWrite<4>(p).last().id = vd2.template getPropRead<2>(q);
+				vd2.template getPropWrite<4>(q).last().id = vd2.template getPropRead<2>(p);
 			}
 
 			++Np;
@@ -1341,6 +1342,10 @@ template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop >
 	vd2.template ghost_put<add_,1>();
 	vd2.template ghost_put<merge_,4>();
 
+#ifdef SE_CLASS3
+	vd2.getDomainIterator();
+#endif
+
 	auto p_it3 = vd.getDomainIterator();
 
 	bool ret = true;
@@ -1348,21 +1353,21 @@ template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop >
 	{
 		auto p = p_it3.get();
 
-		ret &= vd2.template getProp<1>(p) == vd.template getProp<0>(p);
+		ret &= vd2.template getPropRead<1>(p) == vd.template getPropRead<0>(p);
 
 		if (ret == false)
 		{
-			Point<3,float> xp = vd2.getPos(p);
-			std::cout << "ERROR " << vd2.template getProp<1>(p) << "   " << vd.template getProp<0>(p) <<  "    " << xp.toString() << std::endl;
+			Point<3,float> xp = vd2.getPosRead(p);
+			std::cout << "ERROR " << vd2.template getPropWrite<1>(p) << "   " << vd.template getPropWrite<0>(p) <<  "    " << xp.toString() << std::endl;
 		}
 
-		vd.template getProp<3>(p).sort();
-		vd2.template getProp<4>(p).sort();
+		vd.template getPropWrite<3>(p).sort();
+		vd2.template getPropWrite<4>(p).sort();
 
-		ret &= vd.template getProp<3>(p).size() == vd2.template getProp<4>(p).size();
+		ret &= vd.template getPropRead<3>(p).size() == vd2.template getPropRead<4>(p).size();
 
-		for (size_t i = 0 ; i < vd.template getProp<3>(p).size() ; i++)
-			ret &= vd.template getProp<3>(p).get(i).id == vd2.template getProp<4>(p).get(i).id;
+		for (size_t i = 0 ; i < vd.template getPropRead<3>(p).size() ; i++)
+			ret &= vd.template getPropRead<3>(p).get(i).id == vd2.template getPropRead<4>(p).get(i).id;
 
 		if (ret == false)
 			break;
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index c6eb95d385e86442b5426d23827120bfc4aa95de..6c154e19025930bcc4ec187766115a0e6fdd7622 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -1675,13 +1675,13 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 
 			vd.add();
 
-			vd.getLastPos()[0] = key.get(0)*it.getSpacing(0);
-			vd.getLastPos()[1] = key.get(1)*it.getSpacing(1);
-			vd.getLastPos()[2] = key.get(2)*it.getSpacing(2);
+			vd.getLastPosWrite()[0] = key.get(0)*it.getSpacing(0);
+			vd.getLastPosWrite()[1] = key.get(1)*it.getSpacing(1);
+			vd.getLastPosWrite()[2] = key.get(2)*it.getSpacing(2);
 
 			// Fill some properties randomly
 
-			vd.getLastProp<0>() = 0.0;
+			vd.getLastPropWrite<0>() = 0.0;
 
 			++it;
 		}
@@ -1712,12 +1712,12 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 				while (Np.isNext())
 				{
 					auto q = Np.get();
-					Point<3,float> xq = vd.getPos(q);
+					Point<3,float> xq = vd.getPosRead(q);
 
 					float dist = xp.distance(xq);
 
 					if (dist < r_cut)
-						vd.getProp<0>(q) += a*(-dist*dist+r_cut*r_cut);
+						vd.getPropWrite<0>(q) += a*(-dist*dist+r_cut*r_cut);
 
 					++Np;
 				}
@@ -1738,7 +1738,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 				float constant2 = vd.getProp<0>(it3.get());
 				if (fabs(constant - constant2)/constant > eps)
 				{
-					Point<3,float> p = vd.getPos(it3.get());
+					Point<3,float> p = vd.getPosRead(it3.get());
 
 					std::cout << p.toString() << "    " <<  constant2 << "/" << constant << "    " << v_cl.getProcessUnitID() << std::endl;
 					ret = false;
@@ -1755,7 +1755,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 		{
 			auto key = itp.get();
 
-			vd.getProp<0>(key) = 0.0;
+			vd.getPropWrite<0>(key) = 0.0;
 
 			++itp;
 		}
@@ -1772,7 +1772,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 			{
 				// particle p
 				auto p = it2.get();
-				Point<3,float> xp = vd.getPos(p);
+				Point<3,float> xp = vd.getPosRead(p);
 
 				// Get an iterator over the neighborhood particles of p
 				auto Np = NN.getNNIterator<NO_CHECK>(NN.getCell(xp));
@@ -1781,12 +1781,12 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 				while (Np.isNext())
 				{
 					auto q = Np.get();
-					Point<3,float> xq = vd.getPos(q);
+					Point<3,float> xq = vd.getPosRead(q);
 
 					float dist = xp.distance(xq);
 
 					if (dist < r_cut)
-						vd.getProp<0>(q) += a*(-dist*dist+r_cut*r_cut);
+						vd.getPropWrite<0>(q) += a*(-dist*dist+r_cut*r_cut);
 
 					++Np;
 				}
@@ -1799,15 +1799,15 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 			bool ret = true;
 			auto it3 = vd.getDomainIterator();
 
-			float constant = vd.getProp<0>(it3.get());
+			float constant = vd.getPropRead<0>(it3.get());
 			float eps = 0.001;
 
 			while (it3.isNext())
 			{
-				float constant2 = vd.getProp<0>(it3.get());
+				float constant2 = vd.getPropRead<0>(it3.get());
 				if (fabs(constant - constant2)/constant > eps)
 				{
-					Point<3,float> p = vd.getPos(it3.get());
+					Point<3,float> p = vd.getPosRead(it3.get());
 
 					std::cout << p.toString() << "    " <<  constant2 << "/" << constant << "    " << v_cl.getProcessUnitID() << std::endl;
 					ret = false;