diff --git a/regtest/basic/rt-make-0/config b/regtest/basic/rt-make-0/config
index df1f95bf3e..3e5d21664c 100644
--- a/regtest/basic/rt-make-0/config
+++ b/regtest/basic/rt-make-0/config
@@ -1 +1,2 @@
 type=make
+#should we rename this to something meaningful? like "rt-VectorTensor"?
diff --git a/regtest/basic/rt-make-4/main.cpp b/regtest/basic/rt-make-4/main.cpp
index 86f79df350..84315d99a0 100644
--- a/regtest/basic/rt-make-4/main.cpp
+++ b/regtest/basic/rt-make-4/main.cpp
@@ -50,7 +50,9 @@ int main () {
     PLMD::TensorGeneric<4,4> mat;
     PLMD::TensorGeneric<1,4> evec;
     PLMD::VectorGeneric<4> eval_underlying;
-
+    
+    //The both of the two following lines are scary, but on my machine are equivalent
+    //PLMD::Vector1d* eval=reinterpret_cast<PLMD::Vector1d*>(eval_underlying.data());
     auto eval = new(&eval_underlying[0]) PLMD::VectorGeneric<1>;
 
     mat[1][0]=mat[0][1]=3.0;
diff --git a/src/tools/Communicator.h b/src/tools/Communicator.h
index 16ef37ddfb..1a2bf02d1c 100644
--- a/src/tools/Communicator.h
+++ b/src/tools/Communicator.h
@@ -73,13 +73,13 @@ class Communicator {
 /// Init from reference
     template <typename T> explicit Data(T&p): pointer(&p), size(1), nbytes(sizeof(T)), type(getMPIType<T>()) {}
 /// Init from pointer to VectorGeneric
-    template <unsigned n> explicit Data(VectorGeneric<n> *p,int s): pointer(p), size(n*s), nbytes(sizeof(double)), type(getMPIType<double>()) {}
+    template <typename T, unsigned n> explicit Data(VectorTyped<T,n> *p,int s): pointer(p), size(n*s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
 /// Init from reference to VectorGeneric
-    template <unsigned n> explicit Data(VectorGeneric<n> &p): pointer(&p), size(n), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-/// Init from pointer to TensorGeneric
-    template <unsigned n,unsigned m> explicit Data(TensorGeneric<n,m> *p,int s): pointer(p), size(n*m*s), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-/// Init from reference to TensorGeneric
-    template <unsigned n,unsigned m> explicit Data(TensorGeneric<n,m> &p): pointer(&p), size(n*m), nbytes(sizeof(double)), type(getMPIType<double>()) {}
+    template <typename T, unsigned n> explicit Data(VectorTyped<T,n> &p): pointer(&p), size(n), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+/// Init from pointer to TensorTyped
+    template <typename T, unsigned n,unsigned m> explicit Data(TensorTyped<T,n,m> *p,int s): pointer(p), size(n*m*s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+/// Init from reference to TensorTyped
+    template <typename T, unsigned n,unsigned m> explicit Data(TensorTyped<T,n,m> &p): pointer(&p), size(n*m), nbytes(sizeof(T)), type(getMPIType<T>()) {}
 /// Init from reference to std::vector
     template <typename T> explicit Data(std::vector<T>&v) {
       Data d(v.data(),v.size());
@@ -121,10 +121,10 @@ class Communicator {
     MPI_Datatype type;
     template <typename T> explicit ConstData(const T*p,int s): pointer(p), size(s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
     template <typename T> explicit ConstData(const T&p): pointer(&p), size(1), nbytes(sizeof(T)), type(getMPIType<T>()) {}
-    template <unsigned n> explicit ConstData(const VectorGeneric<n> *p,int s): pointer(p), size(n*s), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-    template <unsigned n> explicit ConstData(const VectorGeneric<n> &p): pointer(&p), size(n), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-    template <unsigned n,unsigned m> explicit ConstData(const TensorGeneric<n,m> *p,int s): pointer(p), size(n*m*s), nbytes(sizeof(double)), type(getMPIType<double>()) {}
-    template <unsigned n,unsigned m> explicit ConstData(const TensorGeneric<n,m> &p): pointer(&p), size(n*m), nbytes(sizeof(double)), type(getMPIType<double>()) {}
+    template <typename T,unsigned n> explicit ConstData(const VectorTyped<T,n> *p,int s): pointer(p), size(n*s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+    template <typename T,unsigned n> explicit ConstData(const VectorTyped<T,n> &p): pointer(&p), size(n), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+    template <typename T, unsigned n,unsigned m> explicit ConstData(const TensorTyped<T,n,m> *p,int s): pointer(p), size(n*m*s), nbytes(sizeof(T)), type(getMPIType<T>()) {}
+    template <typename T, unsigned n,unsigned m> explicit ConstData(const TensorTyped<T,n,m> &p): pointer(&p), size(n*m), nbytes(sizeof(T)), type(getMPIType<T>()) {}
     template <typename T> explicit ConstData(const std::vector<T>&v) {
       ConstData d(v.data(),v.size());
       pointer=d.pointer;
diff --git a/src/tools/LoopUnroller.h b/src/tools/LoopUnroller.h
index 6d16761da0..ad1308a3cf 100644
--- a/src/tools/LoopUnroller.h
+++ b/src/tools/LoopUnroller.h
@@ -47,114 +47,89 @@ Implementation is made using template metaprogramming, that is:
 Here xxx is any of the methods of the class.
 
 */
-template<unsigned n>
+//the typename is the second parameter argument so that it can be deduced
+template<typename T, unsigned n>
 class LoopUnroller {
 public:
 /// Set to zero.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]=0.0;`
-  static void _zero(double*d);
+  constexpr static void _zero(T*d);
 /// Add v to d.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]+=v[i];`
-  static void _add(double*d,const double*v);
+  constexpr static void _add(T*d,const T*v);
 /// Subtract v from d.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]-=v[i];`
-  static void _sub(double*d,const double*v);
+  constexpr static void _sub(T*d,const T*v);
 /// Multiply d by s.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]*=s;`
-  static void _mul(double*d,const double s);
+  constexpr static void _mul(T*d,const T s);
 /// Set d to -v.
 /// Same as `for(unsigned i=0;i<n;i++) d[i]=-v[i];`
-  static void _neg(double*d,const double*v);
+  constexpr static void _neg(T*d,const T*v);
 /// Squared modulo of d;
 /// Same as `r=0.0; for(unsigned i=0;i<n;i++) r+=d[i]*d[i]; return r;`
-  static double _sum2(const double*d);
+  constexpr static T _sum2(const T*d);
 /// Dot product of d and v
 /// Same as `r=0.0; for(unsigned i=0;i<n;i++) r+=d[i]*v[i]; return r;`
-  static double _dot(const double*d,const double*v);
+  constexpr static T _dot(const T*d,const T*v);
 };
 
-template<unsigned n>
-void LoopUnroller<n>::_zero(double*d) {
-  LoopUnroller<n-1>::_zero(d);
-  d[n-1]=0.0;
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_zero(T*d) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_zero(d);
+  }
+  d[n-1]=T(0);
 }
 
-template<>
-inline
-void LoopUnroller<1>::_zero(double*d) {
-  d[0]=0.0;
-}
-
-template<unsigned n>
-void LoopUnroller<n>::_add(double*d,const double*a) {
-  LoopUnroller<n-1>::_add(d,a);
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_add(T*d,const T*a) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_add(d,a);
+  }
   d[n-1]+=a[n-1];
 }
 
-template<>
-inline
-void LoopUnroller<1>::_add(double*d,const double*a) {
-  d[0]+=a[0];
-}
-
-template<unsigned n>
-void LoopUnroller<n>::_sub(double*d,const double*a) {
-  LoopUnroller<n-1>::_sub(d,a);
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_sub(T*d,const T*a) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_sub(d,a);
+  }
   d[n-1]-=a[n-1];
 }
 
-template<>
-inline
-void LoopUnroller<1>::_sub(double*d,const double*a) {
-  d[0]-=a[0];
-}
-
-template<unsigned n>
-void LoopUnroller<n>::_mul(double*d,const double s) {
-  LoopUnroller<n-1>::_mul(d,s);
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_mul(T*d,const T s) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_mul(d,s);
+  }
   d[n-1]*=s;
 }
 
-template<>
-inline
-void LoopUnroller<1>::_mul(double*d,const double s) {
-  d[0]*=s;
-}
-
-template<unsigned n>
-void LoopUnroller<n>::_neg(double*d,const double*a ) {
-  LoopUnroller<n-1>::_neg(d,a);
+template<typename T, unsigned n>
+constexpr void LoopUnroller<T,n>::_neg(T*d,const T*a ) {
+  if constexpr (n>1) {
+    LoopUnroller<T,n-1>::_neg(d,a);
+  }
   d[n-1]=-a[n-1];
 }
 
-template<>
-inline
-void LoopUnroller<1>::_neg(double*d,const double*a) {
-  d[0]=-a[0];
+template<typename T, unsigned n>
+constexpr T LoopUnroller<T,n>::_sum2(const T*d) {
+  if constexpr (n>1) {
+    return LoopUnroller<T,n-1>::_sum2(d)+d[n-1]*d[n-1];
+  } else {
+    return d[0]*d[0];
+  }
 }
-
-template<unsigned n>
-double LoopUnroller<n>::_sum2(const double*d) {
-  return LoopUnroller<n-1>::_sum2(d)+d[n-1]*d[n-1];
-}
-
-template<>
-inline
-double LoopUnroller<1>::_sum2(const double*d) {
-  return d[0]*d[0];
-}
-
-template<unsigned n>
-double LoopUnroller<n>::_dot(const double*d,const double*v) {
-  return LoopUnroller<n-1>::_dot(d,v)+d[n-1]*v[n-1];
-}
-
-template<>
-inline
-double LoopUnroller<1>::_dot(const double*d,const double*v) {
-  return d[0]*v[0];
-}
-
+template<typename T, unsigned n>
+constexpr T LoopUnroller<T,n>::_dot(const T*d,const T*v) {
+  if constexpr (n>1) {
+    return LoopUnroller<T,n-1>::_dot(d,v)+d[n-1]*v[n-1];
+  } else {
+    return d[0]*v[0];
+  }
 }
+} //PLMD
 
-#endif
+#endif //__PLUMED_tools_LoopUnroller_h
diff --git a/src/tools/MatrixSquareBracketsAccess.h b/src/tools/MatrixSquareBracketsAccess.h
index 8de3e28f85..1509a8b7ab 100644
--- a/src/tools/MatrixSquareBracketsAccess.h
+++ b/src/tools/MatrixSquareBracketsAccess.h
@@ -73,10 +73,10 @@ class MatrixSquareBracketsAccess {
     // the user should not manipulate it directly
     const MatrixSquareBracketsAccess& t;
     const I i;
-    Const_row(const MatrixSquareBracketsAccess&t,I i); // constructor is private and cannot be manipulated by the user
+    constexpr Const_row(const MatrixSquareBracketsAccess&t,I i); // constructor is private and cannot be manipulated by the user
   public:
     /// access element
-    const C & operator[] (J j)const;
+    constexpr const C & operator[] (J j)const;
   };
 /// Small utility class which just contains a pointer to the T and the row number
   class Row {
@@ -84,28 +84,28 @@ class MatrixSquareBracketsAccess {
     // the user should not manipulate it directly
     MatrixSquareBracketsAccess& t;
     const I i;
-    Row(MatrixSquareBracketsAccess&t,I i); // constructor is private and cannot be manipulated by the user
+    constexpr Row(MatrixSquareBracketsAccess&t,I i); // constructor is private and cannot be manipulated by the user
   public:
     /// access element
-    C & operator[] (J j);
+    constexpr C & operator[] (J j);
   };
 public:
 /// access element (with [][] syntax)
-  Row operator[] (I i);
+  constexpr Row operator[] (I i);
 /// access element (with [][] syntax)
-  Const_row operator[] (I i)const;
+  constexpr Const_row operator[] (I i)const;
 };
 
 template<class T,class C,class I,class J>
-MatrixSquareBracketsAccess<T,C,I,J>::Const_row::Const_row(const MatrixSquareBracketsAccess&t,I i):
+constexpr MatrixSquareBracketsAccess<T,C,I,J>::Const_row::Const_row(const MatrixSquareBracketsAccess&t,I i):
   t(t),i(i) {}
 
 template<class T,class C,class I,class J>
-MatrixSquareBracketsAccess<T,C,I,J>::Row::Row(MatrixSquareBracketsAccess&t,I i):
+constexpr MatrixSquareBracketsAccess<T,C,I,J>::Row::Row(MatrixSquareBracketsAccess&t,I i):
   t(t),i(i) {}
 
 template<class T,class C,class I,class J>
-const C & MatrixSquareBracketsAccess<T,C,I,J>::Const_row::operator[] (J j)const {
+constexpr const C & MatrixSquareBracketsAccess<T,C,I,J>::Const_row::operator[] (J j)const {
 // This appears as a reference to a temporary object
 // however, in reality we know it is a reference to an object that is stored in the
 // t object. We thus suppress the warning raised by cppcheck
@@ -113,7 +113,7 @@ const C & MatrixSquareBracketsAccess<T,C,I,J>::Const_row::operator[] (J j)const
 }
 
 template<class T,class C,class I,class J>
-C & MatrixSquareBracketsAccess<T,C,I,J>::Row::operator[] (J j) {
+constexpr C & MatrixSquareBracketsAccess<T,C,I,J>::Row::operator[] (J j) {
 // This appears as a reference to a temporary object
 // however, in reality we know it is a reference to an object that is stored in the
 // t object. We thus suppress the warning raised by cppcheck
@@ -121,12 +121,12 @@ C & MatrixSquareBracketsAccess<T,C,I,J>::Row::operator[] (J j) {
 }
 
 template<class T,class C,class I,class J>
-typename MatrixSquareBracketsAccess<T,C,I,J>::Row MatrixSquareBracketsAccess<T,C,I,J>::operator[] (I i) {
+constexpr typename MatrixSquareBracketsAccess<T,C,I,J>::Row MatrixSquareBracketsAccess<T,C,I,J>::operator[] (I i) {
   return Row(*this,i);
 }
 
 template<class T,class C,class I,class J>
-typename MatrixSquareBracketsAccess<T,C,I,J>::Const_row MatrixSquareBracketsAccess<T,C,I,J>::operator[] (I i)const {
+constexpr typename MatrixSquareBracketsAccess<T,C,I,J>::Const_row MatrixSquareBracketsAccess<T,C,I,J>::operator[] (I i)const {
   return Const_row(*this,i);
 }
 
diff --git a/src/tools/Tensor.h b/src/tools/Tensor.h
index 4ed456b3d7..6bb813c336 100644
--- a/src/tools/Tensor.h
+++ b/src/tools/Tensor.h
@@ -31,10 +31,10 @@
 
 namespace PLMD {
 
+//should I add something for calling ssyevr?
 /// Small class to contain local utilities.
 /// Should not be used outside of the TensorGeneric class.
-class TensorGenericAux {
-public:
+struct TensorGenericAux {
 // local redefinition, just to avoid including lapack.h here
   static void local_dsyevr(const char *jobz, const char *range, const char *uplo, int *n,
                            double *a, int *lda, double *vl, double *vu, int *
@@ -83,166 +83,146 @@ int main(){
 
 \endverbatim
 */
-template <unsigned n,unsigned m>
-class TensorGeneric:
-  public MatrixSquareBracketsAccess<TensorGeneric<n,m>,double> {
-  std::array<double,n*m> d;
+template<typename T, unsigned n, unsigned m> class TensorTyped;
+
+template<typename T, unsigned n, unsigned m>
+std::ostream & operator<<(std::ostream &os, const TensorTyped<T,n,m>& t);
+
+template<typename T, unsigned n, unsigned m>
+class TensorTyped:
+  public MatrixSquareBracketsAccess<TensorTyped<T,n,m>,T> {
+  std::array<T,n*m> d;
 /// Auxiliary private function for constructor
-  void auxiliaryConstructor();
+  constexpr void auxiliaryConstructor();
 /// Auxiliary private function for constructor
   template<typename... Args>
-  void auxiliaryConstructor(double first,Args... arg);
+  constexpr void auxiliaryConstructor(T first,Args... arg);
+  template <int n_, int m_>
+  using is3x3 = std::enable_if_t<(n_==3 && m_==3)>;
 public:
-/// Constructor accepting n*m double parameters.
+/// Constructor accepting n*m T parameters.
 /// Can be used as Tensor<2,2>(1.0,2.0,3.0,4.0)
 /// In case a wrong number of parameters is given, a static assertion will fail.
   template<typename... Args>
-  TensorGeneric(double first,Args... arg);
+  constexpr TensorTyped(T first,Args... arg);
 /// initialize the tensor to zero
-  TensorGeneric();
+  constexpr TensorTyped();
 /// initialize a tensor as an external product of two Vector
-  TensorGeneric(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2);
+  constexpr TensorTyped(const VectorTyped<T,n>&v1,const VectorTyped<T,m>&v2);
 /// set it to zero
-  void zero();
+  constexpr void zero();
+/// get the underline pointer to data
+  constexpr T* data();
+/// get the underline pointer to data
+  constexpr const T* data() const;
 /// access element
-  double & operator() (unsigned i,unsigned j);
+  constexpr T & operator() (unsigned i,unsigned j);
 /// access element
-  const double & operator() (unsigned i,unsigned j)const;
+  constexpr const T & operator() (unsigned i,unsigned j)const;
 /// increment
-  TensorGeneric& operator +=(const TensorGeneric<n,m>& b);
+  constexpr TensorTyped& operator +=(const TensorTyped& b);
 /// decrement
-  TensorGeneric& operator -=(const TensorGeneric<n,m>& b);
+  constexpr TensorTyped& operator -=(const TensorTyped& b);
 /// multiply
-  TensorGeneric& operator *=(double);
+  constexpr TensorTyped& operator *=(T);
 /// divide
-  TensorGeneric& operator /=(double);
+  constexpr TensorTyped& operator /=(T);
 /// return +t
-  TensorGeneric operator +()const;
+  constexpr TensorTyped operator +()const;
 /// return -t
-  TensorGeneric operator -()const;
+  constexpr TensorTyped operator -()const;
 /// set j-th column
-  TensorGeneric& setCol(unsigned j,const VectorGeneric<n> & c);
+  constexpr TensorTyped& setCol(unsigned j,const VectorTyped<T,n> & c);
 /// set i-th row
-  TensorGeneric& setRow(unsigned i,const VectorGeneric<m> & r);
+  constexpr TensorTyped& setRow(unsigned i,const VectorTyped<T,m> & r);
 /// get j-th column
-  VectorGeneric<n> getCol(unsigned j)const;
+  constexpr VectorTyped<T,n> getCol(unsigned j)const;
 /// get i-th row
-  VectorGeneric<m> getRow(unsigned i)const;
-/// return t1+t2
-  template<unsigned n_,unsigned m_>
-  friend TensorGeneric<n_,m_> operator+(const TensorGeneric<n_,m_>&,const TensorGeneric<n_,m_>&);
-/// return t1+t2
-  template<unsigned n_,unsigned m_>
-  friend TensorGeneric<n_,m_> operator-(const TensorGeneric<n_,m_>&,const TensorGeneric<n_,m_>&);
-/// scale the tensor by a factor s
-  template<unsigned n_,unsigned m_>
-  friend TensorGeneric<n_,m_> operator*(double,const TensorGeneric<n_,m_>&);
-/// scale the tensor by a factor s
-  template<unsigned n_,unsigned m_>
-  friend TensorGeneric<n_,m_> operator*(const TensorGeneric<n_,m_>&,double s);
-/// scale the tensor by a factor 1/s
-  template<unsigned n_,unsigned m_>
-  friend TensorGeneric<n_,m_> operator/(const TensorGeneric<n_,m_>&,double s);
+  constexpr VectorTyped<T,m> getRow(unsigned i)const;
 /// returns the determinant
-  double determinant()const;
+  template<int m_=m, int n_=n,typename = is3x3<m_,n_>>
+  constexpr T determinant()const;
 /// return an identity tensor
-  static TensorGeneric<n,n> identity();
+  static constexpr TensorTyped<T,n,n> identity();
 /// return the matrix inverse
-  TensorGeneric inverse()const;
+  template<int m_=m, int n_=n,typename = is3x3<m_,n_>>
+  constexpr TensorTyped inverse()const;
 /// return the transpose matrix
-  TensorGeneric<m,n> transpose()const;
-/// matrix-matrix multiplication
-  template<unsigned n_,unsigned m_,unsigned l_>
-  friend TensorGeneric<n_,l_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&);
-/// matrix-vector multiplication
-  template<unsigned n_,unsigned m_>
-  friend VectorGeneric<n_> matmul(const TensorGeneric<n_,m_>&,const VectorGeneric<m_>&);
-/// vector-matrix multiplication
-  template<unsigned n_,unsigned m_>
-  friend VectorGeneric<n_> matmul(const VectorGeneric<m_>&,const TensorGeneric<m_,n_>&);
-/// vector-vector multiplication (maps to dotProduct)
-  template<unsigned n_>
-  friend double matmul(const VectorGeneric<n_>&,const VectorGeneric<n_>&);
-/// matrix-matrix-matrix multiplication
-  template<unsigned n_,unsigned m_,unsigned l_,unsigned i_>
-  friend TensorGeneric<n_,i_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&,const TensorGeneric<l_,i_>&);
-/// matrix-matrix-vector multiplication
-  template<unsigned n_,unsigned m_,unsigned l_>
-  friend VectorGeneric<n_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&,const VectorGeneric<l_>&);
-/// vector-matrix-matrix multiplication
-  template<unsigned n_,unsigned m_,unsigned l_>
-  friend VectorGeneric<l_> matmul(const VectorGeneric<n_>&,const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&);
-/// vector-matrix-vector multiplication
-  template<unsigned n_,unsigned m_>
-  friend double matmul(const VectorGeneric<n_>&,const TensorGeneric<n_,m_>&,const VectorGeneric<m_>&);
-/// returns the determinant of a tensor
-  friend double determinant(const TensorGeneric<3,3>&);
-/// returns the inverse of a tensor (same as inverse())
-  friend TensorGeneric<3,3> inverse(const TensorGeneric<3,3>&);
-/// returns the transpose of a tensor (same as transpose())
-  template<unsigned n_,unsigned m_>
-  friend TensorGeneric<n_,m_> transpose(const TensorGeneric<m_,n_>&);
-/// returns the transpose of a tensor (same as TensorGeneric(const VectorGeneric&,const VectorGeneric&))
-  template<unsigned n_,unsigned m_>
-  friend TensorGeneric<n_,m_> extProduct(const VectorGeneric<n>&,const VectorGeneric<m>&);
-  friend TensorGeneric<3,3> dcrossDv1(const VectorGeneric<3>&,const VectorGeneric<3>&);
-  friend TensorGeneric<3,3> dcrossDv2(const VectorGeneric<3>&,const VectorGeneric<3>&);
-  friend TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
-  friend TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&,const VectorGeneric<3>&);
-/// Derivative of a normalized vector
-  friend TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
+  constexpr TensorTyped<T,m,n> transpose()const;
 /// << operator.
 /// Allows printing tensor `t` with `std::cout<<t;`
-  template<unsigned n_,unsigned m_>
-  friend std::ostream & operator<<(std::ostream &os, const TensorGeneric<n_,m_>&);
-/// Diagonalize tensor.
-/// Syntax is the same as Matrix::diagMat.
-/// In addition, it is possible to call if with m_ smaller than n_. In this case,
-/// only the first (smaller) m_ eigenvalues and eigenvectors are retrieved.
-/// If case lapack fails (info!=0) it throws an exception.
-/// Notice that tensor is assumed to be symmetric!!!
-  template<unsigned n_,unsigned m_>
-  friend void diagMatSym(const TensorGeneric<n_,n_>&,VectorGeneric<m_>&evals,TensorGeneric<m_,n_>&evec);
+  friend std::ostream & operator<< <>(std::ostream &os, const TensorTyped&);
 };
 
-template <unsigned n,unsigned m>
-void TensorGeneric<n,m>::auxiliaryConstructor()
+template<typename T, unsigned n, unsigned m>
+constexpr void TensorTyped<T,n,m>::auxiliaryConstructor()
 {}
 
-template <unsigned n,unsigned m>
+template<typename T, unsigned n, unsigned m>
 template<typename... Args>
-void TensorGeneric<n,m>::auxiliaryConstructor(double first,Args... arg) {
+constexpr void TensorTyped<T,n,m>::auxiliaryConstructor(T first,Args... arg) {
   d[n*m-(sizeof...(Args))-1]=first;
   auxiliaryConstructor(arg...);
 }
 
-template <unsigned n,unsigned m>
+template<typename T, unsigned n, unsigned m>
 template<typename... Args>
-TensorGeneric<n,m>::TensorGeneric(double first,Args... arg) {
+constexpr TensorTyped<T,n,m>::TensorTyped(T first,Args... arg) {
   static_assert((sizeof...(Args))+1==n*m,"you are trying to initialize a Tensor with the wrong number of arguments");
   auxiliaryConstructor(first,arg...);
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m>::TensorGeneric() {
-  LoopUnroller<n*m>::_zero(d.data());
+template<typename T, unsigned n, unsigned m>
+constexpr T* TensorTyped<T,n,m>::data() {
+  return d.data();
+}
+template<typename T, unsigned n, unsigned m>
+constexpr const T* TensorTyped<T,n,m>::data() const {
+  return d.data();
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m>::TensorGeneric(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2) {
-  for(unsigned i=0; i<n; i++)
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m>::TensorTyped() {
+  LoopUnroller<T,n*m>::_zero(d.data());
+}
+
+/* between RVO and compile time this should be faster, but slows down openACC, a lot
+template <typename T, unsigned i, unsigned j, unsigned m>
+void external_rec(T*const out,const T*const v1, const T*const v2){
+      if constexpr (j>0) {
+        external_rec<i,j-1,m>(out,v1,v2);
+      } else if constexpr (i>0) {
+          external_rec<i-1,m-1,m>(out,v1,v2);
+      }
+      out[i*m+j]=v1[i]*v2[j];
+}
+template<typename T, unsigned n, unsigned m>
+std::array<T,n*m> externaProd(const VectorGeneric<T,n>&v1,const VectorGeneric<T,m>&v2){
+std::array<T,n*m> toRet;
+external_rec<n-1,m-1,m>(toRet.data(),v1.data(),v2.data());
+return toRet;
+}
+template<typename T, unsigned n, unsigned m>
+TensorGeneric<T,n,m>::TensorGeneric(const VectorGeneric<T,n>&v1,const VectorGeneric<T,m>&v2)
+:d(externaProd(v1,v2)) {}
+*/
+
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m>::TensorTyped(const VectorTyped<T,n>&v1,const VectorTyped<T,m>&v2) {
+  for(unsigned i=0; i<n; i++) {
     for(unsigned j=0; j<m; j++) {
       d[i*m+j]=v1[i]*v2[j];
     }
+  }
 }
 
-template<unsigned n,unsigned m>
-void TensorGeneric<n,m>::zero() {
-  LoopUnroller<n*m>::_zero(d.data());
+template<typename T, unsigned n, unsigned m>
+constexpr void TensorTyped<T,n,m>::zero() {
+  LoopUnroller<T,n*m>::_zero(d.data());
 }
 
-template<unsigned n,unsigned m>
-double & TensorGeneric<n,m>::operator() (unsigned i,unsigned j) {
+template<typename T, unsigned n, unsigned m>
+constexpr T & TensorTyped<T,n,m>::operator() (unsigned i,unsigned j) {
 #ifdef _GLIBCXX_DEBUG
 // index i is implicitly checked by the std::array class
   plumed_assert(j<m);
@@ -250,8 +230,8 @@ double & TensorGeneric<n,m>::operator() (unsigned i,unsigned j) {
   return d[m*i+j];
 }
 
-template<unsigned n,unsigned m>
-const double & TensorGeneric<n,m>::operator() (unsigned i,unsigned j)const {
+template<typename T, unsigned n, unsigned m>
+constexpr const T & TensorTyped<T,n,m>::operator() (unsigned i,unsigned j)const {
 #ifdef _GLIBCXX_DEBUG
 // index i is implicitly checked by the std::array class
   plumed_assert(j<m);
@@ -259,110 +239,110 @@ const double & TensorGeneric<n,m>::operator() (unsigned i,unsigned j)const {
   return d[m*i+j];
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m>& TensorGeneric<n,m>::operator +=(const TensorGeneric<n,m>& b) {
-  LoopUnroller<n*m>::_add(d.data(),b.d.data());
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::operator +=(const TensorTyped<T,n,m>& b) {
+  LoopUnroller<T,n*m>::_add(d.data(),b.d.data());
   return *this;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m>& TensorGeneric<n,m>::operator -=(const TensorGeneric<n,m>& b) {
-  LoopUnroller<n*m>::_sub(d.data(),b.d.data());
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::operator -=(const TensorTyped<T,n,m>& b) {
+  LoopUnroller<T,n*m>::_sub(d.data(),b.d.data());
   return *this;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m>& TensorGeneric<n,m>::operator *=(double s) {
-  LoopUnroller<n*m>::_mul(d.data(),s);
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::operator *=(T s) {
+  LoopUnroller<T,n*m>::_mul(d.data(),s);
   return *this;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m>& TensorGeneric<n,m>::operator /=(double s) {
-  LoopUnroller<n*m>::_mul(d.data(),1.0/s);
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::operator /=(T s) {
+  LoopUnroller<T,n*m>::_mul(d.data(),1.0/s);
   return *this;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m> TensorGeneric<n,m>::operator+()const {
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m> TensorTyped<T,n,m>::operator+()const {
   return *this;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m> TensorGeneric<n,m>::operator-()const {
-  TensorGeneric<n,m> r;
-  LoopUnroller<n*m>::_neg(r.d.data(),d.data());
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m> TensorTyped<T,n,m>::operator-()const {
+  TensorTyped<T,n,m> r;
+  LoopUnroller<T,n*m>::_neg(r.d.data(),d.data());
   return r;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m>& TensorGeneric<n,m>::setCol(unsigned j,const VectorGeneric<n> & c) {
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::setCol(unsigned j,const VectorTyped<T,n> & c) {
   for(unsigned i=0; i<n; ++i) {
     (*this)(i,j)=c(i);
   }
   return *this;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m>& TensorGeneric<n,m>::setRow(unsigned i,const VectorGeneric<m> & r) {
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::setRow(unsigned i,const VectorTyped<T,m> & r) {
   for(unsigned j=0; j<m; ++j) {
     (*this)(i,j)=r(j);
   }
   return *this;
 }
 
-template<unsigned n,unsigned m>
-VectorGeneric<n> TensorGeneric<n,m>::getCol(unsigned j)const {
-  VectorGeneric<n> v;
+template<typename T, unsigned n, unsigned m>
+constexpr VectorTyped<T,n> TensorTyped<T,n,m>::getCol(unsigned j)const {
+  VectorTyped<T,n> v;
   for(unsigned i=0; i<n; ++i) {
     v(i)=(*this)(i,j);
   }
   return v;
 }
 
-template<unsigned n,unsigned m>
-VectorGeneric<m> TensorGeneric<n,m>::getRow(unsigned i)const {
-  VectorGeneric<m> v;
+template<typename T, unsigned n, unsigned m>
+constexpr VectorTyped<T,m> TensorTyped<T,n,m>::getRow(unsigned i)const {
+  VectorTyped<T,m> v;
   for(unsigned j=0; j<m; ++j) {
     v(j)=(*this)(i,j);
   }
   return v;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m> operator+(const TensorGeneric<n,m>&t1,const TensorGeneric<n,m>&t2) {
-  TensorGeneric<n,m> t(t1);
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m> operator+(const TensorTyped<T,n,m>&t1,const TensorTyped<T,n,m>&t2) {
+  TensorTyped<T,n,m> t(t1);
   t+=t2;
   return t;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m> operator-(const TensorGeneric<n,m>&t1,const TensorGeneric<n,m>&t2) {
-  TensorGeneric<n,m> t(t1);
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m> operator-(const TensorTyped<T,n,m>&t1,const TensorTyped<T,n,m>&t2) {
+  TensorTyped<T,n,m> t(t1);
   t-=t2;
   return t;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m> operator*(const TensorGeneric<n,m>&t1,double s) {
-  TensorGeneric<n,m> t(t1);
+template<typename T, typename J, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m> operator*(const TensorTyped<T,n,m>&t1,J s) {
+  TensorTyped<T,n,m> t(t1);
   t*=s;
   return t;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m> operator*(double s,const TensorGeneric<n,m>&t1) {
+template<typename T, typename J, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m> operator*(J s,const TensorTyped<T,n,m>&t1) {
   return t1*s;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m> operator/(const TensorGeneric<n,m>&t1,double s) {
-  return t1*(1.0/s);
+template<typename T, typename J, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m> operator/(const TensorTyped<T,n,m>&t1,J s) {
+  return t1*(T(1.0)/s);
 }
 
-template<>
-inline
-double TensorGeneric<3,3>::determinant()const {
+template<typename T, unsigned n, unsigned m>
+template<int m_, int n_,typename>
+constexpr inline T TensorTyped<T,n,m>::determinant()const {
   return
     d[0]*d[4]*d[8]
     + d[1]*d[5]*d[6]
@@ -372,19 +352,19 @@ double TensorGeneric<3,3>::determinant()const {
     - d[2]*d[4]*d[6];
 }
 
-template<unsigned n,unsigned m>
-inline
-TensorGeneric<n,n> TensorGeneric<n,m>::identity() {
-  TensorGeneric<n,n> t;
+//consider to make this a constexpr function
+template<typename T, unsigned n, unsigned m>
+constexpr  inline TensorTyped<T,n,n> TensorTyped<T,n,m>::identity() {
+  TensorTyped<T,n,n> t;
   for(unsigned i=0; i<n; i++) {
     t(i,i)=1.0;
   }
   return t;
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<m,n> TensorGeneric<n,m>::transpose()const {
-  TensorGeneric<m,n> t;
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,m,n> TensorTyped<T,n,m>::transpose()const {
+  TensorTyped<T,m,n> t;
   for(unsigned i=0; i<m; i++)
     for(unsigned j=0; j<n; j++) {
       t(i,j)=(*this)(j,i);
@@ -392,11 +372,11 @@ TensorGeneric<m,n> TensorGeneric<n,m>::transpose()const {
   return t;
 }
 
-template<>
-inline
-TensorGeneric<3,3> TensorGeneric<3,3>::inverse()const {
-  TensorGeneric t;
-  double invdet=1.0/determinant();
+template<typename T, unsigned n, unsigned m>
+template<int m_, int n_,typename>
+constexpr inline TensorTyped<T,n,m> TensorTyped<T,n,m>::inverse()const {
+  TensorTyped t;
+  T invdet=1.0/determinant();
   for(unsigned i=0; i<3; i++)
     for(unsigned j=0; j<3; j++)
       t(j,i)=invdet*(   (*this)((i+1)%3,(j+1)%3)*(*this)((i+2)%3,(j+2)%3)
@@ -404,9 +384,10 @@ TensorGeneric<3,3> TensorGeneric<3,3>::inverse()const {
   return t;
 }
 
-template<unsigned n,unsigned m,unsigned l>
-TensorGeneric<n,l> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b) {
-  TensorGeneric<n,l> t;
+/// matrix-matrix multiplication
+template<typename T, unsigned n, unsigned m, unsigned l>
+constexpr TensorTyped<T,n,l> matmul(const TensorTyped<T,n,m>&a,const TensorTyped<T,m,l>&b) {
+  TensorTyped<T,n,l> t;
   for(unsigned i=0; i<n; i++)
     for(unsigned j=0; j<l; j++)
       for(unsigned k=0; k<m; k++) {
@@ -415,9 +396,10 @@ TensorGeneric<n,l> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b)
   return t;
 }
 
-template<unsigned n,unsigned m>
-VectorGeneric<n> matmul(const TensorGeneric<n,m>&a,const VectorGeneric<m>&b) {
-  VectorGeneric<n> t;
+/// matrix-vector multiplication
+template<typename T, unsigned n, unsigned m>
+constexpr VectorTyped<T,n> matmul(const TensorTyped<T,n,m>&a,const VectorTyped<T,m>&b) {
+  VectorTyped<T,n> t;
   for(unsigned i=0; i<n; i++)
     for(unsigned j=0; j<m; j++) {
       t(i)+=a(i,j)*b(j);
@@ -425,9 +407,10 @@ VectorGeneric<n> matmul(const TensorGeneric<n,m>&a,const VectorGeneric<m>&b) {
   return t;
 }
 
-template<unsigned n,unsigned m>
-VectorGeneric<n> matmul(const VectorGeneric<m>&a,const TensorGeneric<m,n>&b) {
-  VectorGeneric<n> t;
+/// vector-matrix multiplication
+template<typename T, unsigned n, unsigned m>
+constexpr VectorTyped<T,n> matmul(const VectorTyped<T,m>&a,const TensorTyped<T,m,n>&b) {
+  VectorTyped<T,n> t;
   for(unsigned i=0; i<n; i++)
     for(unsigned j=0; j<m; j++) {
       t(i)+=a(j)*b(j,i);
@@ -435,71 +418,82 @@ VectorGeneric<n> matmul(const VectorGeneric<m>&a,const TensorGeneric<m,n>&b) {
   return t;
 }
 
-template<unsigned n_>
-double matmul(const VectorGeneric<n_>&a,const VectorGeneric<n_>&b) {
+/// vector-vector multiplication (maps to dotProduct)
+template<typename T, unsigned n_>
+constexpr T matmul(const VectorTyped<T,n_>&a,const VectorTyped<T,n_>&b) {
   return dotProduct(a,b);
 }
 
-template<unsigned n,unsigned m,unsigned l,unsigned i>
-TensorGeneric<n,i> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b,const TensorGeneric<l,i>&c) {
+/// matrix-matrix-matrix multiplication
+template<typename T, unsigned n, unsigned m, unsigned l, unsigned i>
+constexpr TensorTyped<T,n,i> matmul(const TensorTyped<T,n,m>&a,const TensorTyped<T,m,l>&b,const TensorTyped<T,l,i>&c) {
   return matmul(matmul(a,b),c);
 }
 
-template<unsigned n,unsigned m,unsigned l>
-VectorGeneric<n> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b,const VectorGeneric<l>&c) {
+/// matrix-matrix-vector multiplication
+template<typename T, unsigned n, unsigned m, unsigned l>
+constexpr VectorTyped<T,n> matmul(const TensorTyped<T,n,m>&a,const TensorTyped<T,m,l>&b,const VectorTyped<T,l>&c) {
   return matmul(matmul(a,b),c);
 }
 
-template<unsigned n,unsigned m,unsigned l>
-VectorGeneric<l> matmul(const VectorGeneric<n>&a,const TensorGeneric<n,m>&b,const TensorGeneric<m,l>&c) {
+/// vector-matrix-matrix multiplication
+template<typename T, unsigned n, unsigned m, unsigned l>
+constexpr VectorTyped<T,l> matmul(const VectorTyped<T,n>&a,const TensorTyped<T,n,m>&b,const TensorTyped<T,m,l>&c) {
   return matmul(matmul(a,b),c);
 }
 
-template<unsigned n,unsigned m>
-double matmul(const VectorGeneric<n>&a,const TensorGeneric<n,m>&b,const VectorGeneric<m>&c) {
+/// vector-matrix-vector multiplication
+template<typename T, unsigned n, unsigned m>
+constexpr T matmul(const VectorTyped<T,n>&a,const TensorTyped<T,n,m>&b,const VectorTyped<T,m>&c) {
   return matmul(matmul(a,b),c);
 }
 
+template <typename T>
 inline
-double determinant(const TensorGeneric<3,3>&t) {
+constexpr T determinant(const TensorTyped<T,3,3>&t) {
   return t.determinant();
 }
 
+template <typename T>
 inline
-TensorGeneric<3,3> inverse(const TensorGeneric<3,3>&t) {
+constexpr TensorTyped<T,3,3> inverse(const TensorTyped<T,3,3>&t) {
   return t.inverse();
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m> transpose(const TensorGeneric<m,n>&t) {
+/// returns the transpose of a tensor (same as TensorGeneric::transpose())
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m> transpose(const TensorTyped<T,m,n>&t) {
   return t.transpose();
 }
 
-template<unsigned n,unsigned m>
-TensorGeneric<n,m> extProduct(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2) {
-  return TensorGeneric<n,m>(v1,v2);
+/// returns the transpose of a tensor (same as TensorGeneric(const VectorGeneric&,const VectorGeneric&))
+template<typename T, unsigned n, unsigned m>
+constexpr TensorTyped<T,n,m> extProduct(const VectorTyped<T,n>&v1,const VectorTyped<T,m>&v2) {
+  return TensorTyped<T,n,m>(v1,v2);
 }
 
+template <typename T>
 inline
-TensorGeneric<3,3> dcrossDv1(const VectorGeneric<3>&v1,const VectorGeneric<3>&v2) {
+constexpr TensorTyped<T,3,3> dcrossDv1(const VectorTyped<T,3>&v1,const VectorTyped<T,3>&v2) {
   (void) v1; // this is to avoid warnings. still the syntax of this function is a bit dummy...
-  return TensorGeneric<3,3>(
-           0.0, v2[2],-v2[1],
-           -v2[2],   0.0, v2[0],
-           v2[1],-v2[0],   0.0);
+  return TensorTyped<T,3,3>(
+           T(0.0),  v2[2], -v2[1],
+           -v2[2], T(0.0), v2[0],
+           v2[1],  -v2[0], T(0.0));
 }
 
+template <typename T>
 inline
-TensorGeneric<3,3> dcrossDv2(const VectorGeneric<3>&v1,const VectorGeneric<3>&v2) {
+constexpr TensorTyped<T,3,3> dcrossDv2(const VectorTyped<T,3>&v1,const VectorTyped<T,3>&v2) {
   (void) v2; // this is to avoid warnings. still the syntax of this function is a bit dummy...
-  return TensorGeneric<3,3>(
-           0.0,-v1[2],v1[1],
-           v1[2],0.0,-v1[0],
-           -v1[1],v1[0],0.0);
+  return TensorTyped<T,3,3>(
+           T(0.0),-v1[2], v1[1],
+           v1[2], T(0.0),-v1[0],
+           -v1[1],v1[0], T(0.0));
 }
 
-template<unsigned n,unsigned m>
-std::ostream & operator<<(std::ostream &os, const TensorGeneric<n,m>& t) {
+template<typename T, unsigned n, unsigned m>
+std::ostream & operator<<(std::ostream &os, const TensorTyped<T,n,m>& t) {
   for(unsigned i=0; i<n; i++)
     for(unsigned j=0; j<m; j++) {
       if(i!=(n-1) || j!=(m-1)) {
@@ -511,6 +505,9 @@ std::ostream & operator<<(std::ostream &os, const TensorGeneric<n,m>& t) {
   return os;
 }
 
+template <unsigned n,unsigned m>
+using TensorGeneric=TensorTyped<double,n,m>;
+
 /// \ingroup TOOLBOX
 typedef TensorGeneric<1,1> Tensor1d;
 /// \ingroup TOOLBOX
@@ -524,33 +521,40 @@ typedef TensorGeneric<5,5> Tensor5d;
 /// \ingroup TOOLBOX
 typedef Tensor3d Tensor;
 
+template <typename T>
 inline
-TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2) {
-
-  TensorGeneric<3,3> t;
+TensorTyped<T,3,3> VcrossTensor(const VectorTyped<T,3>&v1,const TensorTyped<T,3,3>&v2) {
+  TensorTyped<T,3,3> t;
   for(unsigned i=0; i<3; i++) {
     t.setRow(i,matmul(dcrossDv2(v1,v1),v2.getRow(i)));
   }
   return t;
 }
 
+template <typename T>
 inline
-TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&v2,const VectorGeneric<3>&v1) {
-  TensorGeneric<3,3> t;
+TensorTyped<T,3,3> VcrossTensor(const TensorTyped<T,3,3>&v2,const VectorTyped<T,3>&v1) {
+  TensorTyped<T,3,3> t;
   for(unsigned i=0; i<3; i++) {
     t.setRow(i,-matmul(dcrossDv2(v1,v1),v2.getRow(i)));
   }
   return t;
 }
 
-
+template <typename T>
 inline
-TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2) {
+TensorTyped<T,3,3> deriNorm(const VectorTyped<T,3>&v1,const TensorTyped<T,3,3>&v2) {
   // delta(v) = delta(v1/v1.norm) = 1/v1.norm*(delta(v1) - (v.delta(v1))cross v;
-  double over_norm = 1./v1.modulo();
+  T over_norm = 1./v1.modulo();
   return over_norm*(v2 - over_norm*over_norm*(extProduct(matmul(v2,v1),v1)));
 }
 
+/// Diagonalize tensor.
+/// Syntax is the same as Matrix::diagMat.
+/// In addition, it is possible to call if with m_ smaller than n_. In this case,
+/// only the first (smaller) m_ eigenvalues and eigenvectors are retrieved.
+/// If case lapack fails (info!=0) it throws an exception.
+/// Notice that tensor is assumed to be symmetric!!!
 template<unsigned n,unsigned m>
 void diagMatSym(const TensorGeneric<n,n>&mat,VectorGeneric<m>&evals,TensorGeneric<m,n>&evec) {
   // some guess number to make sure work is large enough.
@@ -573,8 +577,8 @@ void diagMatSym(const TensorGeneric<n,n>&mat,VectorGeneric<m>&evals,TensorGeneri
   int info=0;            // result
   int liwork=iwork.size();
   int lwork=work.size();
-  TensorGenericAux::local_dsyevr("V", (n==m?"A":"I"), "U", &nn, const_cast<double*>(&mat_copy[0][0]), &nn, &vl, &vu, &one, &mm,
-                                 &abstol, &mout, &evals_tmp[0], &evec[0][0], &nn,
+  TensorGenericAux::local_dsyevr("V", (n==m?"A":"I"), "U", &nn, mat_copy.data(), &nn, &vl, &vu, &one, &mm,
+                                 &abstol, &mout, evals_tmp.data(), evec.data(), &nn,
                                  isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
   if(info!=0)
     plumed_error()<<"Error diagonalizing matrix\n"
@@ -604,7 +608,6 @@ void diagMatSym(const TensorGeneric<n,n>&mat,VectorGeneric<m>&evals,TensorGeneri
 
 static_assert(sizeof(Tensor)==9*sizeof(double), "code cannot work if this is not satisfied");
 
-}
+} //PLMD
 
 #endif
-
diff --git a/src/tools/Tools.h b/src/tools/Tools.h
index a1d07795f9..3142984084 100644
--- a/src/tools/Tools.h
+++ b/src/tools/Tools.h
@@ -85,7 +85,7 @@ class Tools {
   static bool convertToInt(const std::string & str,T &t);
 /// @brief the  recursive part of the template fastpow implementation
   template <int exp, typename T=double, std::enable_if_t< (exp >=0), bool> = true>
-  static inline /*consteval*/ T fastpow_rec(T base, T result);
+  static constexpr inline T fastpow_rec(T base, T result);
 public:
 /// Split the line in words using separators.
 /// It also take into account parenthesis. Outer parenthesis found are removed from
@@ -178,10 +178,11 @@ class Tools {
 /// extension("pippo/.t")="" whereas extension("pippo/a.t")="t"
   static std::string extension(const std::string&);
 /// Fast int power
-  static double fastpow(double base,int exp);
+  template <typename T>
+  static constexpr inline T fastpow(T base,int exp);
 /// Fast int power for power known at compile time
   template <int exp, typename T=double>
-  static inline /*consteval*/ T fastpow(T base);
+  static constexpr inline T fastpow(T base);
 /// Modified 0th-order Bessel function of the first kind
   static double bessel0(const double& val);
 /// Check if a string full starts with string start.
@@ -248,16 +249,16 @@ class Tools {
     if(s==0) {
       return;
     }
-    set_to_zero(&vec[0][0],s*n);
+    set_to_zero(vec[0].data(),s*n);
   }
 
-  template<unsigned n,unsigned m>
-  static void set_to_zero(std::vector<TensorGeneric<n,m>> & vec) {
+  template<typename T, unsigned n,unsigned m>
+  static void set_to_zero(std::vector<TensorTyped<T,n,m>> & vec) {
     unsigned s=vec.size();
     if(s==0) {
       return;
     }
-    set_to_zero(&vec[0](0,0),s*n*m);
+    set_to_zero(vec[0].data(),s*n*m);
   }
 
   static std::unique_ptr<std::lock_guard<std::mutex>> molfile_lock();
@@ -487,8 +488,8 @@ bool Tools::convertNoexcept(T i,std::string & str) {
   return true;
 }
 
-inline
-double Tools::fastpow(double base, int exp) {
+template <typename T>
+constexpr inline T Tools::fastpow(T base, int exp) {
   if(exp<0) {
     exp=-exp;
     base=1.0/base;
@@ -506,7 +507,7 @@ double Tools::fastpow(double base, int exp) {
 }
 
 template <int exp, typename T, std::enable_if_t< (exp >=0), bool>>
-inline T Tools::fastpow_rec(T const base, T result) {
+constexpr inline T Tools::fastpow_rec(T const base, T result) {
   if constexpr (exp == 0) {
     return result;
   }
@@ -517,7 +518,7 @@ inline T Tools::fastpow_rec(T const base, T result) {
 }
 
 template <int exp, typename T>
-inline T Tools::fastpow(T const base) {
+constexpr inline T Tools::fastpow(T const base) {
   if constexpr (exp<0) {
     return  fastpow_rec<-exp,T>(1.0/base,1.0);
   } else {
diff --git a/src/tools/Vector.h b/src/tools/Vector.h
index 0f585be328..97fd3b59bd 100644
--- a/src/tools/Vector.h
+++ b/src/tools/Vector.h
@@ -73,239 +73,260 @@ int main(){
 \endverbatim
 
 */
+template<typename T, unsigned n> class VectorTyped;
 
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n> delta(const VectorTyped<T, n>&,const VectorTyped<T, n>&);
 
-template <unsigned n>
-class VectorGeneric {
-  std::array<double,n> d;
+template<typename T, unsigned n>
+constexpr T dotProduct(const VectorTyped<T, n>&,const VectorTyped<T, n>&);
+
+
+template<typename T, unsigned n>
+std::ostream & operator<<(std::ostream &os, const VectorTyped<T, n>& v);
+template<typename T, unsigned n>
+class VectorTyped {
+  std::array<T,n> d;
 /// Auxiliary private function for constructor
-  void auxiliaryConstructor();
+  constexpr void auxiliaryConstructor();
 /// Auxiliary private function for constructor
   template<typename... Args>
-  void auxiliaryConstructor(double first,Args... arg);
+  constexpr void auxiliaryConstructor(T first,Args... arg);
 public:
-/// Constructor accepting n double parameters.
+/// Constructor accepting n T parameters.
 /// Can be used as Vector<3>(1.0,2.0,3.0) or Vector<2>(2.0,3.0).
 /// In case a wrong number of parameters is given, a static assertion will fail.
   template<typename... Args>
-  VectorGeneric(double first,Args... arg);
+  constexpr VectorTyped(T first,Args... arg);
 /// create it null
-  VectorGeneric();
+  constexpr VectorTyped();
+/// get the underline pointer to the data
+  constexpr T* data();
+/// get the underline pointer to the data
+  constexpr const T* data() const;
 /// set it to zero
-  void zero();
+  constexpr void zero();
 /// array-like access [i]
-  double & operator[](unsigned i);
+  constexpr   T & operator[](unsigned i);
 /// array-like access [i]
-  const double & operator[](unsigned i)const;
+  constexpr const T & operator[](unsigned i)const;
 /// parenthesis access (i)
-  double & operator()(unsigned i);
+  constexpr   T & operator()(unsigned i);
 /// parenthesis access (i)
-  const double & operator()(unsigned i)const;
+  constexpr const T & operator()(unsigned i)const;
 /// increment
-  VectorGeneric& operator +=(const VectorGeneric& b);
+  constexpr VectorTyped& operator +=(const VectorTyped& b);
 /// decrement
-  VectorGeneric& operator -=(const VectorGeneric& b);
+  constexpr VectorTyped& operator -=(const VectorTyped& b);
 /// multiply
-  VectorGeneric& operator *=(double s);
+  constexpr VectorTyped& operator *=(T s);
 /// divide
-  VectorGeneric& operator /=(double s);
+  constexpr VectorTyped& operator /=(T s);
 /// sign +
-  VectorGeneric operator +()const;
+  constexpr VectorTyped operator +()const;
 /// sign -
-  VectorGeneric operator -()const;
+  constexpr VectorTyped operator -()const;
 /// return v1+v2
-  template<unsigned m>
-  friend VectorGeneric<m> operator+(const VectorGeneric<m>&,const VectorGeneric<m>&);
+  template<typename U, unsigned m>
+  friend constexpr VectorTyped<U, m> operator+(const VectorTyped<U, m>&,const VectorTyped<U, m>&);
 /// return v1-v2
-  template<unsigned m>
-  friend VectorGeneric<m> operator-(const VectorGeneric<m>&,const VectorGeneric<m>&);
+  template<typename U, unsigned m>
+  friend constexpr VectorTyped<U, m> operator-(VectorTyped<U, m>,const VectorTyped<U, m>&);
 /// return s*v
-  template<unsigned m>
-  friend VectorGeneric<m> operator*(double,const VectorGeneric<m>&);
+  template<typename U, typename J, unsigned m>
+  friend constexpr VectorTyped<U, m> operator*(J,VectorTyped<U, m>);
 /// return v*s
-  template<unsigned m>
-  friend VectorGeneric<m> operator*(const VectorGeneric<m>&,double);
+  template<typename U, typename J, unsigned m>
+  friend constexpr VectorTyped<U, m> operator*(VectorTyped<U, m>,J);
 /// return v/s
-  template<unsigned m>
-  friend VectorGeneric<m> operator/(const VectorGeneric<m>&,double);
+  template<typename U, typename J, unsigned m>
+  friend constexpr VectorTyped<U, m> operator/(const VectorTyped<U, m>&,J);
 /// return v2-v1
-  template<unsigned m>
-  friend VectorGeneric<m> delta(const VectorGeneric<m>&v1,const VectorGeneric<m>&v2);
+  friend constexpr VectorTyped delta<>(const VectorTyped&v1,const VectorTyped&v2);
 /// return v1 .scalar. v2
-  template<unsigned m>
-  friend double dotProduct(const VectorGeneric<m>&,const VectorGeneric<m>&);
+  friend constexpr T dotProduct<>(const VectorTyped&,const VectorTyped&);
+  //this bad boy produces a warning (in fact becasue declrare the crossproduc as a friend for ALL thhe possible combinations of n and T)
 /// return v1 .vector. v2
 /// Only available for size 3
-  friend VectorGeneric<3> crossProduct(const VectorGeneric<3>&,const VectorGeneric<3>&);
+  template<typename U>
+  friend constexpr VectorTyped<U, 3> crossProduct(const VectorTyped<U, 3>&,const VectorTyped<U, 3>&);
 /// compute the squared modulo
-  double modulo2()const;
+  constexpr T modulo2()const;
 /// Compute the modulo.
 /// Shortcut for sqrt(v.modulo2())
-  double modulo()const;
+  constexpr T modulo()const;
 /// friend version of modulo2 (to simplify some syntax)
-  template<unsigned m>
-  friend double modulo2(const VectorGeneric<m>&);
+  template<typename U, unsigned m>
+  friend constexpr U modulo2(const VectorTyped<U, m>&);
 /// friend version of modulo (to simplify some syntax)
-  template<unsigned m>
-  friend double modulo(const VectorGeneric<m>&);
+  template<typename U, unsigned m>
+  friend constexpr U modulo(const VectorTyped<U, m>&);
 /// << operator.
 /// Allows printing vector `v` with `std::cout<<v;`
-  template<unsigned m>
-  friend std::ostream & operator<<(std::ostream &os, const VectorGeneric<m>&);
+  friend std::ostream & operator<< <>(std::ostream &os, const VectorTyped&);
 };
 
-template <unsigned n>
-void VectorGeneric<n>::auxiliaryConstructor()
+template<typename T, unsigned n>
+constexpr void VectorTyped<T, n>::auxiliaryConstructor()
 {}
 
-template <unsigned n>
+template<typename T, unsigned n>
 template<typename... Args>
-void VectorGeneric<n>::auxiliaryConstructor(double first,Args... arg) {
+constexpr void VectorTyped<T, n>::auxiliaryConstructor(T first,Args... arg) {
   d[n-(sizeof...(Args))-1]=first;
   auxiliaryConstructor(arg...);
 }
 
-template <unsigned n>
+template<typename T, unsigned n>
 template<typename... Args>
-VectorGeneric<n>::VectorGeneric(double first,Args... arg) {
+constexpr VectorTyped<T, n>::VectorTyped(T first,Args... arg) {
   static_assert((sizeof...(Args))+1==n,"you are trying to initialize a Vector with the wrong number of arguments");
   auxiliaryConstructor(first,arg...);
 }
 
-template <unsigned n>
-void VectorGeneric<n>::zero() {
-  LoopUnroller<n>::_zero(d.data());
+template<typename T, unsigned n>
+constexpr T* VectorTyped<T, n>::data() {
+  return d.data();
+}
+
+template<typename T, unsigned n>
+constexpr const T* VectorTyped<T, n>::data() const {
+  return d.data();
 }
 
-template <unsigned n>
-VectorGeneric<n>::VectorGeneric() {
-  LoopUnroller<n>::_zero(d.data());
+template<typename T, unsigned n>
+constexpr void VectorTyped<T, n>::zero() {
+  LoopUnroller<T, n>::_zero(d.data());
 }
 
-template <unsigned n>
-double & VectorGeneric<n>::operator[](unsigned i) {
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n>::VectorTyped() {
+  LoopUnroller<T, n>::_zero(d.data());
+}
+
+template<typename T, unsigned n>
+constexpr T & VectorTyped<T, n>::operator[](unsigned i) {
   return d[i];
 }
 
-template <unsigned n>
-const double & VectorGeneric<n>::operator[](unsigned i)const {
+template<typename T, unsigned n>
+constexpr const T & VectorTyped<T, n>::operator[](unsigned i)const {
   return d[i];
 }
 
-template <unsigned n>
-double & VectorGeneric<n>::operator()(unsigned i) {
+template<typename T, unsigned n>
+constexpr T & VectorTyped<T, n>::operator()(unsigned i) {
   return d[i];
 }
 
-template <unsigned n>
-const double & VectorGeneric<n>::operator()(unsigned i)const {
+template<typename T, unsigned n>
+constexpr const T & VectorTyped<T, n>::operator()(unsigned i)const {
   return d[i];
 }
 
-template <unsigned n>
-VectorGeneric<n>& VectorGeneric<n>::operator +=(const VectorGeneric<n>& b) {
-  LoopUnroller<n>::_add(d.data(),b.d.data());
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n>& VectorTyped<T, n>::operator +=(const VectorTyped<T, n>& b) {
+  LoopUnroller<T, n>::_add(d.data(),b.d.data());
   return *this;
 }
 
-template <unsigned n>
-VectorGeneric<n>& VectorGeneric<n>::operator -=(const VectorGeneric<n>& b) {
-  LoopUnroller<n>::_sub(d.data(),b.d.data());
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n>& VectorTyped<T, n>::operator -=(const VectorTyped<T, n>& b) {
+  LoopUnroller<T, n>::_sub(d.data(),b.d.data());
   return *this;
 }
 
-template <unsigned n>
-VectorGeneric<n>& VectorGeneric<n>::operator *=(double s) {
-  LoopUnroller<n>::_mul(d.data(),s);
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n>& VectorTyped<T, n>::operator *=(T s) {
+  LoopUnroller<T, n>::_mul(d.data(),s);
   return *this;
 }
 
-template <unsigned n>
-VectorGeneric<n>& VectorGeneric<n>::operator /=(double s) {
-  LoopUnroller<n>::_mul(d.data(),1.0/s);
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n>& VectorTyped<T, n>::operator /=(T s) {
+  LoopUnroller<T, n>::_mul(d.data(),1.0/s);
   return *this;
 }
 
-template <unsigned n>
-VectorGeneric<n>  VectorGeneric<n>::operator +()const {
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n>  VectorTyped<T, n>::operator +()const {
   return *this;
 }
 
-template <unsigned n>
-VectorGeneric<n> VectorGeneric<n>::operator -()const {
-  VectorGeneric<n> r;
-  LoopUnroller<n>::_neg(r.d.data(),d.data());
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n> VectorTyped<T, n>::operator -()const {
+  VectorTyped<T, n> r;
+  LoopUnroller<T, n>::_neg(r.d.data(),d.data());
   return r;
 }
 
-template <unsigned n>
-VectorGeneric<n> operator+(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
-  VectorGeneric<n> v(v1);
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n> operator+(const VectorTyped<T, n>&v1,const VectorTyped<T, n>&v2) {
+  VectorTyped<T, n> v(v1);
   return v+=v2;
 }
 
-template <unsigned n>
-VectorGeneric<n> operator-(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
-  VectorGeneric<n> v(v1);
-  return v-=v2;
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n> operator-(VectorTyped<T, n>v1,const VectorTyped<T, n>&v2) {
+  return v1-=v2;
 }
 
-template <unsigned n>
-VectorGeneric<n> operator*(double s,const VectorGeneric<n>&v) {
-  VectorGeneric<n> vv(v);
-  return vv*=s;
+template<typename T, typename J, unsigned n>
+constexpr VectorTyped<T, n> operator*(J s,VectorTyped<T, n>v) {
+  return v*=s;
 }
 
-template <unsigned n>
-VectorGeneric<n> operator*(const VectorGeneric<n>&v,double s) {
-  return s*v;
+template<typename T, typename J, unsigned n>
+constexpr VectorTyped<T, n> operator*(VectorTyped<T, n> v,J s) {
+  return v*=s;
 }
 
-template <unsigned n>
-VectorGeneric<n> operator/(const VectorGeneric<n>&v,double s) {
-  return v*(1.0/s);
+template<typename T, typename J, unsigned n>
+constexpr VectorTyped<T, n> operator/(const VectorTyped<T, n>&v,J s) {
+  return v*(T(1.0)/s);
 }
 
-template <unsigned n>
-VectorGeneric<n> delta(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
+template<typename T, unsigned n>
+constexpr VectorTyped<T, n> delta(const VectorTyped<T, n>&v1,const VectorTyped<T, n>&v2) {
   return v2-v1;
 }
 
-template <unsigned n>
-double VectorGeneric<n>::modulo2()const {
-  return LoopUnroller<n>::_sum2(d.data());
+template<typename T, unsigned n>
+constexpr T VectorTyped<T, n>::modulo2()const {
+  return LoopUnroller<T, n>::_sum2(d.data());
 }
 
-template <unsigned n>
-double dotProduct(const VectorGeneric<n>& v1,const VectorGeneric<n>& v2) {
-  return LoopUnroller<n>::_dot(v1.d.data(),v2.d.data());
+template<typename T, unsigned n>
+constexpr T dotProduct(const VectorTyped<T, n>& v1,const VectorTyped<T, n>& v2) {
+  return LoopUnroller<T, n>::_dot(v1.d.data(),v2.d.data());
 }
 
-inline
-VectorGeneric<3> crossProduct(const VectorGeneric<3>& v1,const VectorGeneric<3>& v2) {
-  return VectorGeneric<3>(
+template<typename T>
+constexpr inline
+VectorTyped<T, 3> crossProduct(const VectorTyped<T, 3>& v1,const VectorTyped<T, 3>& v2) {
+  return VectorTyped<T, 3>(
            v1[1]*v2[2]-v1[2]*v2[1],
            v1[2]*v2[0]-v1[0]*v2[2],
            v1[0]*v2[1]-v1[1]*v2[0]);
 }
 
-template<unsigned n>
-double VectorGeneric<n>::modulo()const {
+template<typename T, unsigned n>
+constexpr T VectorTyped<T, n>::modulo()const {
   return sqrt(modulo2());
 }
 
-template<unsigned n>
-double modulo2(const VectorGeneric<n>&v) {
+template<typename T, unsigned n>
+constexpr T modulo2(const VectorTyped<T, n>&v) {
   return v.modulo2();
 }
 
-template<unsigned n>
-double modulo(const VectorGeneric<n>&v) {
+template<typename T, unsigned n>
+constexpr T modulo(const VectorTyped<T, n>&v) {
   return v.modulo();
 }
 
-template<unsigned n>
-std::ostream & operator<<(std::ostream &os, const VectorGeneric<n>& v) {
+template<typename T, unsigned n>
+std::ostream & operator<<(std::ostream &os, const VectorTyped<T, n>& v) {
   for(unsigned i=0; i<n-1; i++) {
     os<<v(i)<<" ";
   }
@@ -313,6 +334,8 @@ std::ostream & operator<<(std::ostream &os, const VectorGeneric<n>& v) {
   return os;
 }
 
+template<unsigned n>
+using VectorGeneric=VectorTyped<double, n>;
 
 /// \ingroup TOOLBOX
 /// Alias for one dimensional vectors
@@ -331,13 +354,13 @@ typedef VectorGeneric<4> Vector4d;
 typedef VectorGeneric<5> Vector5d;
 /// \ingroup TOOLBOX
 /// Alias for three dimensional vectors
-typedef Vector3d Vector;
+using Vector=Vector3d;
+//using the using keyword seems to be more trasparent
 
 static_assert(sizeof(VectorGeneric<2>)==2*sizeof(double), "code cannot work if this is not satisfied");
 static_assert(sizeof(VectorGeneric<3>)==3*sizeof(double), "code cannot work if this is not satisfied");
 static_assert(sizeof(VectorGeneric<4>)==4*sizeof(double), "code cannot work if this is not satisfied");
 
-}
-
-#endif
+} //PLMD
 
+#endif //__PLUMED_tools_Vector_h