diff --git a/CHANGELOG.md b/CHANGELOG.md
index f5879995733915ec38daeab5d0cb7ab8c6969807..610fba5d92942a91fed614f5fc9b2bce168d932d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,11 +1,10 @@
 # Change Log
 All notable changes to this project will be documented in this file.
 
-## [0.5.0] - 
+## [0.5.0] - mid July 2016
 
 ### Added
-- map_list map communicate particles across processors mooving the information of all the particle map_list give the possibility to
-           give a list of property to move from one to another processor
+- map_list map communicate particles across processors mooving the information of all the particle map_list give the possibility to give a list of property to move from one to another processor
 - Numeric: Finite Differences discretization with matrix contruction and parallel solvers
 
 ### Fixed
@@ -82,3 +81,18 @@ All notable changes to this project will be documented in this file.
 ### Changed
 - Nothing to report
 
+
+
+# Planned in the next Releases
+
+## [0.7.0] - Mid of October
+
+### Added
+- Dynamic Load Balancies examples and interface fixation
+- Check Point restart
+
+## [0.6.0] - Beginning of september
+
+### Added
+- Parallel IO, new formats, imroved writers
+
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index a3778f022c1ac1be9229ac52f5eb2896c9e7c9f3..0ad8579643db74ed1e98912ab9c26856db8580c4 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -27,6 +27,7 @@
 #include "Vector/vector_dist_ofb.hpp"
 #include "Decomposition/CartDecomposition.hpp"
 #include "data_type/aggregate.hpp"
+#include "vector_dist_operators.hpp"
 
 #define V_SUB_UNIT_FACTOR 64
 
@@ -883,6 +884,21 @@ public:
 		return v_prp.template get<id>(vec_key.getKey());
 	}
 
+	/*! \brief Get the property of an element
+	 *
+	 * see the vector_dist iterator usage to get an element key
+	 *
+	 * \tparam id property id
+	 * \param vec_key vector element
+	 *
+	 * \return return the selected property of the vector element
+	 *
+	 */
+	template<unsigned int id> inline const auto getProp(vect_dist_key_dx vec_key) const -> decltype(v_prp.template get<id>(vec_key.getKey()))
+	{
+		return v_prp.template get<id>(vec_key.getKey());
+	}
+
 	/*! \brief It move all the particles that does not belong to the local processor to the respective processor
 	 *
 	 * \tparam out of bound policy it specify what to do when the particles are detected out of bound
@@ -1438,7 +1454,7 @@ public:
 	 * \return an iterator
 	 *
 	 */
-	vector_dist_iterator getGhostIterator()
+	vector_dist_iterator getGhostIterator() const
 	{
 		return vector_dist_iterator(g_m, v_pos.size());
 	}
@@ -1448,7 +1464,7 @@ public:
 	 * \return an iterator
 	 *
 	 */
-	vector_dist_iterator getDomainIterator()
+	vector_dist_iterator getDomainIterator() const
 	{
 		return vector_dist_iterator(0, g_m);
 	}
@@ -1611,6 +1627,16 @@ public:
 #endif
 		return v_cl;
 	}
+
+	//////////////////////////// Vector expression implementation ////////////////////////
+
+	template <unsigned int prp> inline vector_dist_expression<prp,vector_dist<dim,St,prop,Decomposition,Memory> > getV()
+	{
+		vector_dist_expression<prp,vector_dist<dim,St,prop,Decomposition,Memory> > exp_v(*this);
+
+		return exp_v;
+	}
 };
 
+
 #endif /* VECTOR_HPP_ */
diff --git a/src/Vector/vector_dist_operators.hpp b/src/Vector/vector_dist_operators.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..da74fadac3c5cc6d0625b65a6a9c804248ffb0e3
--- /dev/null
+++ b/src/Vector/vector_dist_operators.hpp
@@ -0,0 +1,650 @@
+/*
+ * vector_dist_operators.hpp
+ *
+ *  Created on: Jun 11, 2016
+ *      Author: i-bird
+ */
+
+#ifndef SRC_VECTOR_VECTOR_DIST_OPERATORS_HPP_
+#define SRC_VECTOR_VECTOR_DIST_OPERATORS_HPP_
+
+#define SUM 1
+#define SUB 2
+#define MUL 3
+#define DIV 4
+
+template <typename exp1, typename exp2, unsigned int op>
+class vector_dist_expression_op
+{
+
+};
+
+template <typename exp1, typename exp2>
+class vector_dist_expression_op<exp1,exp2,SUM>
+{
+	const exp1 o1;
+	const exp2 o2;
+
+public:
+
+	typedef decltype(o1.value(vect_dist_key_dx(0))) r_type;
+
+	vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+	:o1(o1),o2(o2)
+	{}
+
+	template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+	{
+		return o1.value(key) + o2.value(key);
+	}
+};
+
+template <typename exp1, typename exp2>
+class vector_dist_expression_op<exp1,exp2,SUB>
+{
+	const exp1 o1;
+	const exp2 o2;
+
+public:
+
+	vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+	:o1(o1),o2(o2)
+	{}
+
+	template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+	{
+		return o1.value(key) - o2.value(key);
+	}
+
+};
+
+template <typename exp1, typename exp2>
+class vector_dist_expression_op<exp1,exp2,MUL>
+{
+	const exp1 o1;
+	const exp2 o2;
+
+public:
+
+	vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+	:o1(o1),o2(o2)
+	{}
+
+	template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+	{
+		return o1.value(key) * o2.value(key);
+	}
+
+};
+
+template <typename exp1, typename exp2>
+class vector_dist_expression_op<exp1,exp2,DIV>
+{
+	const exp1 o1;
+	const exp2 o2;
+
+public:
+
+	vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
+	:o1(o1),o2(o2)
+	{}
+
+	template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type > inline r_type value(const vect_dist_key_dx & key) const
+	{
+		return o1.value(key) / o2.value(key);
+	}
+
+};
+
+/*! \brief Main class that encapsulate a vector operation expression
+ *
+ * \param prp property involved
+ * \param exp expression type
+ *
+ */
+template<unsigned int prp, typename vector>
+class vector_dist_expression
+{
+	vector & v;
+
+public:
+
+	vector_dist_expression(vector & v)
+	:v(v)
+	{}
+
+	inline auto value(const vect_dist_key_dx & k) const -> decltype(v.template getProp<prp>(k))
+	{
+		return v.template getProp<prp>(k);
+	}
+
+	/*! \brief evaluate the expression and store it in the vector
+	 *
+	 * \param v_exp expression to evaluate
+	 *
+	 */
+	template<typename exp1, typename exp2, unsigned int op> vector & operator=(const vector_dist_expression_op<exp1,exp2,op> & v_exp)
+	{
+		auto it = v.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			v.template getProp<prp>(key) = v_exp.value(key);
+
+			++it;
+		}
+
+		return v;
+	}
+
+	/*! \brief evaluate the expression and store it in the vector
+	 *
+	 * \param v_exp expression to evaluate
+	 *
+	 */
+	vector & operator=(double d)
+	{
+		auto it = v.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			v.template getProp<prp>(key) = d;
+
+			++it;
+		}
+
+		return v;
+	}
+};
+
+/*! \brief Main class that encapsulate a vector operation expression
+ *
+ * \param prp property involved
+ * \param exp expression type
+ *
+ */
+template<unsigned int prp>
+class vector_dist_expression<prp,double>
+{
+	double d;
+
+public:
+
+	inline vector_dist_expression(double & d)
+	:d(d)
+	{}
+
+	inline double value(const vect_dist_key_dx & k) const
+	{
+		return d;
+	}
+};
+
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int p1, unsigned int p2, typename v1, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,SUM>
+operator+(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,SUM> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,SUM>
+operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<prp1,v1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,SUM> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,SUM>
+operator+(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,SUM> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,SUM>
+operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,SUM> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int prp1 , typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,SUM>
+operator+(const vector_dist_expression<prp1,v1> & va, double d)
+{
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,SUM> exp_sum(va,vector_dist_expression<0,double>(d));
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int prp1 , typename v1>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,SUM>
+operator+(double d, const vector_dist_expression<prp1,v1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,SUM> exp_sum(vector_dist_expression<0,double>(d),vb);
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1 , typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,SUM>
+operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,SUM> exp_sum(va,vector_dist_expression<0,double>(d));
+
+	return exp_sum;
+}
+
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int p1, unsigned int p2, typename v1, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,SUB>
+operator-(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,SUB> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1, unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p2,v2>,SUB>
+operator-(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<p2,v2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p2,v2>,SUB> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1, unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression_op<exp1,exp2,op1>,SUB>
+operator-(const vector_dist_expression<p2,v2> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<p2,v2>, vector_dist_expression_op<exp1,exp2,op1>,SUB> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1, typename exp3, typename exp4, unsigned int op2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,SUB>
+operator-(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,SUB> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,SUB>
+operator-(const vector_dist_expression<prp1,v1> & va, double d)
+{
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,SUB> exp_sum(va,vector_dist_expression<0,double>(d));
+
+	return exp_sum;
+}
+
+/* \brief sum two distributed vectors
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,SUB>
+operator-(double d, const vector_dist_expression<prp1,v1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,SUB> exp_sum(vector_dist_expression<0,double>(d),vb);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<p2,v2>,MUL>
+operator*(double d, const vector_dist_expression<p2,v2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<p2,v2>,MUL> exp_sum(vector_dist_expression<0,double>(d),vb);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression<0,double>,MUL>
+operator*(const vector_dist_expression<p2,v2> & va, double d)
+{
+	vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression<0,double>,MUL> exp_sum(va,vector_dist_expression<0,double>(d));
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int p1, typename v1,unsigned int p2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,MUL>
+operator*(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,MUL> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int p1, typename v1, typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression_op<exp1,exp2,op1>,MUL>
+operator*(const vector_dist_expression<p1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression_op<exp1,exp2,op1>,MUL> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int p1, typename v1, typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p1,v1>,MUL>
+operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<p1,v1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p1,v1>,MUL> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,MUL>
+operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,MUL> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,DIV>
+operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,DIV> exp_sum(va,vector_dist_expression<0,double>(d));
+
+	return exp_sum;
+}
+
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1, typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,DIV>
+operator/(double d, const vector_dist_expression_op<exp1,exp2,op1> & va)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,DIV> exp_sum(vector_dist_expression<0,double>(d),va);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,DIV>
+operator/(const vector_dist_expression<prp1,v1> & va, double d)
+{
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,DIV> exp_sum(va,vector_dist_expression<0,double>(d));
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,DIV>
+operator/(double d, const vector_dist_expression<prp1,v1> & va)
+{
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,DIV> exp_sum(va,vector_dist_expression<0,double>(d));
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1, unsigned int prp2, typename v2>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<prp2,v2>,DIV>
+operator/(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression<prp2,v2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<prp2,v2>,DIV> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1, typename exp1,typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,DIV>
+operator/(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,DIV> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<unsigned int prp1, typename v1, typename exp1,typename exp2, unsigned int op1>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,DIV>
+operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<prp1,v1> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,DIV> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+/* \brief multiplication double with a vector expression
+ *
+ * \param v1 vector1
+ * \param v2 vector2
+ *
+ * \return a grid_key_dx_expression that encapsulate the expression
+ *
+ */
+template<typename exp1,typename exp2, unsigned int op1, typename exp3, typename exp4, unsigned int op2>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,DIV>
+operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,DIV> exp_sum(va,vb);
+
+	return exp_sum;
+}
+
+#endif /* SRC_VECTOR_VECTOR_DIST_OPERATORS_HPP_ */
diff --git a/src/Vector/vector_dist_operators_unit_tests.hpp b/src/Vector/vector_dist_operators_unit_tests.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fbbed13fb76f80d862c30777e7a735452132a560
--- /dev/null
+++ b/src/Vector/vector_dist_operators_unit_tests.hpp
@@ -0,0 +1,437 @@
+/*
+ * vector_dist_operators_unit_tests.hpp
+ *
+ *  Created on: Jun 11, 2016
+ *      Author: i-bird
+ */
+
+#ifndef SRC_VECTOR_VECTOR_DIST_OPERATORS_UNIT_TESTS_HPP_
+#define SRC_VECTOR_VECTOR_DIST_OPERATORS_UNIT_TESTS_HPP_
+
+constexpr int A = 0;
+constexpr int B = 1;
+constexpr int C = 2;
+
+constexpr int PA = 3;
+constexpr int PB = 4;
+constexpr int PC = 5;
+
+//////////////////// Here we define all the function to checl the operators
+
+template <unsigned int prp, typename vector> bool check_values(const vector & v,float a)
+{
+	bool ret = true;
+	auto it = v.getDomainIterator();
+
+	while (it.isNext())
+	{
+		ret &= v.template getProp<prp>(it.get()) == a;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_complex_expr(const vector & vd)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd.template getProp<B>(key) + 2.0 + vd.template getProp<B>(key) - 2.0*vd.template getProp<C>(key) / 5.0;
+		float base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_sum(const vector & vd, double d)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd.template getProp<B>(key) + d;
+		float base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_sum(const vector & vd1, const vector & vd2)
+{
+	bool ret = true;
+	auto it = vd1.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd1.template getProp<B>(key) + vd2.template getProp<C>(key);
+		float base2 = vd1.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_sum_3(const vector & vd1)
+{
+	bool ret = true;
+	auto it = vd1.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd1.template getProp<B>(key) + vd1.template getProp<C>(key) + vd1.template getProp<B>(key);
+		float base2 = vd1.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_sum_4(const vector & vd1)
+{
+	bool ret = true;
+	auto it = vd1.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd1.template getProp<B>(key) + vd1.template getProp<C>(key) + vd1.template getProp<B>(key) + vd1.template getProp<C>(key);
+		float base2 = vd1.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_sub(const vector & vd, double d)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd.template getProp<B>(key) - d;
+		float base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_sub(double d, const vector & vd)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = d - vd.template getProp<B>(key);
+		float base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_sub(const vector & vd1, const vector & vd2)
+{
+	bool ret = true;
+	auto it = vd1.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd1.template getProp<C>(key) - vd2.template getProp<B>(key);
+		float base2 = vd1.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_mul(const vector & vd, double d)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd.template getProp<B>(key) * d;
+		float base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+
+template <typename vector> bool check_values_mul(const vector & vd1, const vector & vd2)
+{
+	bool ret = true;
+	auto it = vd1.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd1.template getProp<C>(key) * vd2.template getProp<B>(key);
+		float base2 = vd1.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+
+template <typename vector> bool check_values_div(const vector & vd, double d)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd.template getProp<B>(key) / d;
+		float base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_div(double d, const vector & vd)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = d / vd.template getProp<B>(key);
+		float base2 = vd.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> bool check_values_div(const vector & vd1, const vector & vd2)
+{
+	bool ret = true;
+	auto it = vd1.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		float base1 = vd1.template getProp<C>(key) / vd2.template getProp<B>(key);
+		float base2 = vd1.template getProp<A>(key);
+
+		ret &=  base1 == base2;
+
+		++it;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
+
+template <typename vector> void fill_values(const vector & v)
+{
+	auto it = v.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto p = it.get();
+
+		v.template getProp<A>(p) = p.getKey()+1;
+		v.template getProp<B>(p) = 2.0*p.getKey()+1;
+		v.template getProp<C>(p) = 3.0*p.getKey()+1;
+
+		++it;
+	}
+}
+
+float getProp0(vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>>> & v, size_t p)
+{
+	return v.template getProp<A>(p);
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+BOOST_AUTO_TEST_SUITE( vector_dist_test )
+
+
+BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
+{
+	if (create_vcluster().getProcessingUnits() > 3)
+		return;
+
+	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+	// Boundary conditions
+	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+	// ghost
+	Ghost<3,float> ghost(0.05);
+
+	vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>>> vd(100,box,bc,ghost);
+
+	vd.getV<A>() = 1.0;
+	vd.getV<B>() = 2.0f;
+	vd.getV<C>() = 3.0;
+
+	check_values<A>(vd,1.0);
+	check_values<B>(vd,2.0);
+	check_values<C>(vd,3.0);
+
+	fill_values(vd);
+
+	vd.getV<A>() = vd.getV<B>() + 2.0 + vd.getV<B>() - 2.0*vd.getV<C>() / 5.0;
+	check_values_complex_expr(vd);
+
+	// Various combination of 2 operator
+
+	vd.getV<A>() = vd.getV<B>() + 2.0;
+	check_values_sum(vd,2.0);
+	vd.getV<A>() = 2.0 + vd.getV<B>();
+	check_values_sum(vd,2.0);
+	vd.getV<A>() = vd.getV<C>() + vd.getV<B>();
+	check_values_sum(vd,vd);
+
+	vd.getV<A>() = vd.getV<B>() - 2.0;
+	check_values_sub(vd,2.0);
+	vd.getV<A>() = 2.0 - vd.getV<B>();
+	check_values_sub(2.0,vd);
+	vd.getV<A>() = vd.getV<C>() - vd.getV<B>();
+	check_values_sub(vd,vd);
+
+	vd.getV<A>() = vd.getV<B>() * 2.0;
+	check_values_mul(vd,2.0);
+	vd.getV<A>() = 2.0 * vd.getV<B>();
+	check_values_mul(vd,2.0);
+	vd.getV<A>() = vd.getV<C>() * vd.getV<B>();
+	check_values_mul(vd,vd);
+
+	vd.getV<A>() = vd.getV<B>() / 2.0;
+	check_values_div(vd,2.0);
+	vd.getV<A>() = 2.0 / vd.getV<B>();
+	check_values_div(vd,2.0);
+	vd.getV<A>() = vd.getV<C>() / vd.getV<B>();
+	check_values_div(vd,vd);
+
+	// Variuos combination 3 operator
+
+	vd.getV<A>() = vd.getV<B>() + (vd.getV<C>() + vd.getV<B>());
+	check_values_sum_3(vd);
+	vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) + vd.getV<B>();
+	check_values_sum_3(vd);
+	vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) + (vd.getV<C>() + vd.getV<B>());
+	check_values_sum_4(vd);
+
+/*	vd.getV<A>() = vd.getV<B>() - (vd.getV<C>() + vd.getV<B>());
+	vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) - vd.getV<B>();
+	vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) - (vd.getV<C>() + vd.getV<B>());
+
+	vd.getV<A>() = vd.getV<B>() * (vd.getV<C>() + vd.getV<B>());
+	vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) * vd.getV<B>();
+	vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) * (vd.getV<C>() + vd.getV<B>());
+
+	vd.getV<A>() = vd.getV<B>() / (vd.getV<C>() + vd.getV<B>());
+	vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) / vd.getV<B>();
+	vd.getV<A>() = (vd.getV<C>() + vd.getV<B>()) / (vd.getV<C>() + vd.getV<B>());*/
+
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+
+#endif /* SRC_VECTOR_VECTOR_DIST_OPERATORS_UNIT_TESTS_HPP_ */
diff --git a/src/main.cpp b/src/main.cpp
index 3a554a6e44ba254cfb243b548d73b408f899f7c7..98435d5c9e8c1d69ecb0807b99721d853106b34b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -26,6 +26,7 @@
 #include "dec_optimizer_unit_test.hpp"
 #include "Grid/grid_dist_id_unit_test.hpp"
 #include "Vector/vector_dist_unit_test.hpp"
+#include "Vector/vector_dist_operators_unit_tests.hpp"
 #include "Decomposition/Distribution/Distribution_unit_tests.hpp"
 //#include "DLB/DLB_unit_test.hpp"
 #include "Graph/dist_map_graph_unit_test.hpp"