diff --git a/Grid/simd/Grid_vector_types.h b/Grid/simd/Grid_vector_types.h
index 4f952bb2..fd71a84a 100644
--- a/Grid/simd/Grid_vector_types.h
+++ b/Grid/simd/Grid_vector_types.h
@@ -1055,7 +1055,7 @@ accelerator_inline typename toComplexMapper<Rsimd>::Complexified toComplex(const
 }
 
 
-accelerator_inline void precisionChange(vRealF    *out,vRealD    *in,int nvec)
+accelerator_inline void precisionChange(vRealF    *out,const vRealD    *in,int nvec)
 {
   assert((nvec&0x1)==0);
   for(int m=0;m*2<nvec;m++){
@@ -1063,7 +1063,7 @@ accelerator_inline void precisionChange(vRealF    *out,vRealD    *in,int nvec)
     out[m].v=Optimization::PrecisionChange::DtoS(in[n].v,in[n+1].v);
   }
 }
-accelerator_inline void precisionChange(vRealH    *out,vRealD    *in,int nvec)
+accelerator_inline void precisionChange(vRealH    *out,const vRealD    *in,int nvec)
 {
   assert((nvec&0x3)==0);
   for(int m=0;m*4<nvec;m++){
@@ -1071,7 +1071,7 @@ accelerator_inline void precisionChange(vRealH    *out,vRealD    *in,int nvec)
     out[m].v=Optimization::PrecisionChange::DtoH(in[n].v,in[n+1].v,in[n+2].v,in[n+3].v);
   }
 }
-accelerator_inline void precisionChange(vRealH    *out,vRealF    *in,int nvec)
+accelerator_inline void precisionChange(vRealH    *out,const vRealF    *in,int nvec)
 {
   assert((nvec&0x1)==0);
   for(int m=0;m*2<nvec;m++){
@@ -1079,7 +1079,7 @@ accelerator_inline void precisionChange(vRealH    *out,vRealF    *in,int nvec)
     out[m].v=Optimization::PrecisionChange::StoH(in[n].v,in[n+1].v);
   }
 }
-accelerator_inline void precisionChange(vRealD    *out,vRealF    *in,int nvec)
+accelerator_inline void precisionChange(vRealD    *out,const vRealF    *in,int nvec)
 {
   assert((nvec&0x1)==0);
   for(int m=0;m*2<nvec;m++){
@@ -1095,7 +1095,7 @@ accelerator_inline void precisionChange(vRealD    *out,vRealF    *in,int nvec)
     //  |                                                 ~~~~~~~^
   }
 }
-accelerator_inline void precisionChange(vRealD    *out,vRealH    *in,int nvec)
+accelerator_inline void precisionChange(vRealD    *out,const vRealH    *in,int nvec)
 {
   assert((nvec&0x3)==0);
   for(int m=0;m*4<nvec;m++){
@@ -1103,7 +1103,7 @@ accelerator_inline void precisionChange(vRealD    *out,vRealH    *in,int nvec)
     Optimization::PrecisionChange::HtoD(in[m].v,out[n].v,out[n+1].v,out[n+2].v,out[n+3].v);
   }
 }
-accelerator_inline void precisionChange(vRealF    *out,vRealH    *in,int nvec)
+accelerator_inline void precisionChange(vRealF    *out,const vRealH    *in,int nvec)
 {
   assert((nvec&0x1)==0);
   for(int m=0;m*2<nvec;m++){
@@ -1111,12 +1111,12 @@ accelerator_inline void precisionChange(vRealF    *out,vRealH    *in,int nvec)
     Optimization::PrecisionChange::HtoS(in[m].v,out[n].v,out[n+1].v);
   }
 }
-accelerator_inline void precisionChange(vComplexF *out,vComplexD *in,int nvec){ precisionChange((vRealF *)out,(vRealD *)in,nvec);}
-accelerator_inline void precisionChange(vComplexH *out,vComplexD *in,int nvec){ precisionChange((vRealH *)out,(vRealD *)in,nvec);}
-accelerator_inline void precisionChange(vComplexH *out,vComplexF *in,int nvec){ precisionChange((vRealH *)out,(vRealF *)in,nvec);}
-accelerator_inline void precisionChange(vComplexD *out,vComplexF *in,int nvec){ precisionChange((vRealD *)out,(vRealF *)in,nvec);}
-accelerator_inline void precisionChange(vComplexD *out,vComplexH *in,int nvec){ precisionChange((vRealD *)out,(vRealH *)in,nvec);}
-accelerator_inline void precisionChange(vComplexF *out,vComplexH *in,int nvec){ precisionChange((vRealF *)out,(vRealH *)in,nvec);}
+accelerator_inline void precisionChange(vComplexF *out,const vComplexD *in,int nvec){ precisionChange((vRealF *)out,(vRealD *)in,nvec);}
+accelerator_inline void precisionChange(vComplexH *out,const vComplexD *in,int nvec){ precisionChange((vRealH *)out,(vRealD *)in,nvec);}
+accelerator_inline void precisionChange(vComplexH *out,const vComplexF *in,int nvec){ precisionChange((vRealH *)out,(vRealF *)in,nvec);}
+accelerator_inline void precisionChange(vComplexD *out,const vComplexF *in,int nvec){ precisionChange((vRealD *)out,(vRealF *)in,nvec);}
+accelerator_inline void precisionChange(vComplexD *out,const vComplexH *in,int nvec){ precisionChange((vRealD *)out,(vRealH *)in,nvec);}
+accelerator_inline void precisionChange(vComplexF *out,const vComplexH *in,int nvec){ precisionChange((vRealF *)out,(vRealH *)in,nvec);}
 
 // Check our vector types are of an appropriate size.
 
diff --git a/Grid/simd/Grid_vector_unops.h b/Grid/simd/Grid_vector_unops.h
index b89bb785..35f1721c 100644
--- a/Grid/simd/Grid_vector_unops.h
+++ b/Grid/simd/Grid_vector_unops.h
@@ -29,8 +29,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 			   /*  END LEGAL */
-#ifndef GRID_VECTOR_UNOPS
-#define GRID_VECTOR_UNOPS
+#pragma once
 
 #include <cmath>
 
@@ -112,6 +111,9 @@ template <class scalar>
 struct ImagFunctor {
   accelerator scalar operator()(const scalar &a) const { return imag(a); }
 };
+/////////////
+// Unary operations
+/////////////
 template <class S, class V>
 accelerator_inline Grid_simd<S, V> real(const Grid_simd<S, V> &r) {
   return SimdApply(RealFunctor<S>(), r);
@@ -168,6 +170,65 @@ template <class S, class V>
 accelerator_inline Grid_simd<S, V> div(const Grid_simd<S, V> &r, Integer y) {
   return SimdApply(DivIntFunctor<S>(y), r);
 }
+/// Double 2 cases
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> real(const Grid_simd2<S, V> &r) {
+  return SimdApply(RealFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> imag(const Grid_simd2<S, V> &r) {
+  return SimdApply(ImagFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> sqrt(const Grid_simd2<S, V> &r) {
+  return SimdApply(SqrtRealFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> cos(const Grid_simd2<S, V> &r) {
+  return SimdApply(CosRealFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> sin(const Grid_simd2<S, V> &r) {
+  return SimdApply(SinRealFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> acos(const Grid_simd2<S, V> &r) {
+  return SimdApply(AcosRealFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> asin(const Grid_simd2<S, V> &r) {
+  return SimdApply(AsinRealFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> log(const Grid_simd2<S, V> &r) {
+  return SimdApply(LogRealFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> abs(const Grid_simd2<S, V> &r) {
+  return SimdApply(AbsRealFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> exp(const Grid_simd2<S, V> &r) {
+  return SimdApply(ExpFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> Not(const Grid_simd2<S, V> &r) {
+  return SimdApply(NotFunctor<S>(), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> pow(const Grid_simd2<S, V> &r, double y) {
+  return SimdApply(PowRealFunctor<S>(y), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> mod(const Grid_simd2<S, V> &r, Integer y) {
+  return SimdApply(ModIntFunctor<S>(y), r);
+}
+template <class S, class V>
+accelerator_inline Grid_simd2<S, V> div(const Grid_simd2<S, V> &r, Integer y) {
+  return SimdApply(DivIntFunctor<S>(y), r);
+}
+
+
 ////////////////////////////////////////////////////////////////////////////
 // Allows us to assign into **conformable** real vectors from complex
 ////////////////////////////////////////////////////////////////////////////
@@ -193,23 +254,22 @@ struct OrOrFunctor {
 ////////////////////////////////
 template <class S, class V>
 accelerator_inline Grid_simd<S, V> operator&(const Grid_simd<S, V> &x,
-				 const Grid_simd<S, V> &y) {
+					     const Grid_simd<S, V> &y) {
   return SimdApplyBinop(AndFunctor<S>(), x, y);
 }
 template <class S, class V>
 accelerator_inline Grid_simd<S, V> operator&&(const Grid_simd<S, V> &x,
-				  const Grid_simd<S, V> &y) {
+					      const Grid_simd<S, V> &y) {
   return SimdApplyBinop(AndAndFunctor<S>(), x, y);
 }
 template <class S, class V>
 accelerator_inline Grid_simd<S, V> operator|(const Grid_simd<S, V> &x,
-				 const Grid_simd<S, V> &y) {
+					     const Grid_simd<S, V> &y) {
   return SimdApplyBinop(OrFunctor<S>(), x, y);
 }
 template <class S, class V>
 accelerator_inline Grid_simd<S, V> operator||(const Grid_simd<S, V> &x,
-				  const Grid_simd<S, V> &y) {
+					      const Grid_simd<S, V> &y) {
   return SimdApplyBinop(OrOrFunctor<S>(), x, y);
 }
 NAMESPACE_END(Grid);
-#endif
diff --git a/Grid/simd/Simd.h b/Grid/simd/Simd.h
index 76ca3bef..ddee6a36 100644
--- a/Grid/simd/Simd.h
+++ b/Grid/simd/Simd.h
@@ -224,18 +224,14 @@ accelerator_inline void Gpermute(VectorSIMD &y,const VectorSIMD &b,int perm);
 NAMESPACE_END(Grid);
 
 #include <Grid/simd/Grid_vector_types.h>
+#include <Grid/simd/Grid_doubled_vector.h>
 #include <Grid/simd/Grid_vector_unops.h>
 
 NAMESPACE_BEGIN(Grid);
-// Default precision
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
+
+// Default precision is wired to double
 typedef vRealD vReal;
 typedef vComplexD vComplex;
-#else
-typedef vRealF vReal;
-typedef vComplexF vComplex;
-#endif
-
  
 inline std::ostream& operator<< (std::ostream& stream, const vComplexF &o){
   int nn=vComplexF::Nsimd();