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 875bb4f5f0..47ae0673f7 100644
--- a/src/tools/Communicator.h
+++ b/src/tools/Communicator.h
@@ -76,10 +76,10 @@ class Communicator {
     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 <typename T, unsigned n> explicit Data(VectorTyped<T,n> &p): pointer(&p), size(n), nbytes(sizeof(T)), type(getMPIType<T>()) {}
-/// 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>()) {}
+/// 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()); pointer=d.pointer; size=d.size; type=d.type;
@@ -106,8 +106,8 @@ class Communicator {
     template <typename T> explicit ConstData(const T&p): pointer(&p), size(1), nbytes(sizeof(T)), type(getMPIType<T>()) {}
     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 <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,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; size=d.size; type=d.type;
     }
diff --git a/src/tools/Tensor.h b/src/tools/Tensor.h
index c831bc7607..755050204e 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,142 @@ int main(){
 
 \endverbatim
 */
-template <unsigned n,unsigned m>
-class TensorGeneric:
-  public MatrixSquareBracketsAccess<TensorGeneric<n,m>,double>
+
+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<double,n*m> d;
+  std::array<T,n*m> d;
 /// Auxiliary private function for constructor
   void auxiliaryConstructor();
 /// Auxiliary private function for constructor
   template<typename... Args>
-  void auxiliaryConstructor(double first,Args... arg);
+  void auxiliaryConstructor(T first,Args... arg);
 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);
+  TensorTyped(T first,Args... arg);
 /// initialize the tensor to zero
-  TensorGeneric();
+  TensorTyped();
 /// initialize a tensor as an external product of two Vector
-  TensorGeneric(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2);
+  TensorTyped(const VectorTyped<T,n>&v1,const VectorTyped<T,m>&v2);
 /// set it to zero
   void zero();
+/// get the underline pointer to data
+  T* data();
+/// get the underline pointer to data
+  const T* data() const;
 /// access element
-  double & operator() (unsigned i,unsigned j);
+  T & operator() (unsigned i,unsigned j);
 /// access element
-  const double & operator() (unsigned i,unsigned j)const;
+  const T & operator() (unsigned i,unsigned j)const;
 /// increment
-  TensorGeneric& operator +=(const TensorGeneric<n,m>& b);
+  TensorTyped& operator +=(const TensorTyped& b);
 /// decrement
-  TensorGeneric& operator -=(const TensorGeneric<n,m>& b);
+  TensorTyped& operator -=(const TensorTyped& b);
 /// multiply
-  TensorGeneric& operator *=(double);
+  TensorTyped& operator *=(T);
 /// divide
-  TensorGeneric& operator /=(double);
+  TensorTyped& operator /=(T);
 /// return +t
-  TensorGeneric operator +()const;
+  TensorTyped operator +()const;
 /// return -t
-  TensorGeneric operator -()const;
+  TensorTyped operator -()const;
 /// set j-th column
-  TensorGeneric& setCol(unsigned j,const VectorGeneric<n> & c);
+  TensorTyped& setCol(unsigned j,const VectorTyped<T,n> & c);
 /// set i-th row
-  TensorGeneric& setRow(unsigned i,const VectorGeneric<m> & r);
+  TensorTyped& setRow(unsigned i,const VectorTyped<T,m> & r);
 /// get j-th column
-  VectorGeneric<n> getCol(unsigned j)const;
+  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);
+  VectorTyped<T,m> getRow(unsigned i)const;
 /// returns the determinant
-  double determinant()const;
+  T determinant()const;
 /// return an identity tensor
-  static TensorGeneric<n,n> identity();
+  static TensorTyped<T,n,n> identity();
 /// return the matrix inverse
-  TensorGeneric inverse()const;
+  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>&);
+  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>
+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)
+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)
+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<double,n*m>::_zero(d.data());
+template<typename T, unsigned n, unsigned m>
+T* TensorTyped<T,n,m>::data() {return d.data();}
+template<typename T, unsigned n, unsigned m>
+const T* TensorTyped<T,n,m>::data() const {return d.data();}
+
+template<typename T, unsigned n, unsigned m>
+TensorTyped<T,n,m>::TensorTyped() {
+  LoopUnroller<T,n*m>::_zero(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++)for(unsigned j=0; j<m; j++)d[i*m+j]=v1[i]*v2[j];
+/* 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<unsigned n,unsigned m>
-void TensorGeneric<n,m>::zero() {
-  LoopUnroller<double,n*m>::_zero(d.data());
+template<typename T, unsigned n, unsigned m>
+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>
-double & TensorGeneric<n,m>::operator() (unsigned i,unsigned j) {
+template<typename T, unsigned n, unsigned m>
+void TensorTyped<T,n,m>::zero() {
+  LoopUnroller<T,n*m>::_zero(d.data());
+}
+
+template<typename T, unsigned n, unsigned m>
+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 +226,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>
+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,102 +235,103 @@ 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<double,n*m>::_add(d.data(),b.d.data());
+template<typename T, unsigned n, unsigned m>
+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<double,n*m>::_sub(d.data(),b.d.data());
+template<typename T, unsigned n, unsigned m>
+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<double,n*m>::_mul(d.data(),s);
+template<typename T, unsigned n, unsigned m>
+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<double,n*m>::_mul(d.data(),1.0/s);
+template<typename T, unsigned n, unsigned m>
+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>
+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<double,n*m>::_neg(r.d.data(),d.data());
+template<typename T, unsigned n, unsigned m>
+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>
+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>
+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>
+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>
+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>
+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>
+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>
+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>
+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>
+TensorTyped<T,n,m> operator/(const TensorTyped<T,n,m>&t1,J s) {
+  return t1*(T(1.0)/s);
 }
 
-template<>
+template<typename T, unsigned n, unsigned m>
 inline
-double TensorGeneric<3,3>::determinant()const {
+T TensorTyped<T,n,m>::determinant()const {
+  static_assert(n==3&&m==3,"determinanat can be called only for 3x3 Tensors");
   return
     d[0]*d[4]*d[8]
     + d[1]*d[5]*d[6]
@@ -364,120 +341,136 @@ double TensorGeneric<3,3>::determinant()const {
     - d[2]*d[4]*d[6];
 }
 
-template<unsigned n,unsigned m>
+//consider to make this a constexpr function
+template<typename T, unsigned n, unsigned m>
 inline
-TensorGeneric<n,n> TensorGeneric<n,m>::identity() {
-  TensorGeneric<n,n> t;
+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>
+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);
   return t;
 }
 
-template<>
+template<typename T, unsigned n, unsigned m>
 inline
-TensorGeneric<3,3> TensorGeneric<3,3>::inverse()const {
-  TensorGeneric t;
-  double invdet=1.0/determinant();
+TensorTyped<T,n,m> TensorTyped<T,n,m>::inverse()const {
+  static_assert(n==3&&m==3,"inverse can be called only for 3x3 Tensors");
+  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)
                         -(*this)((i+1)%3,(j+2)%3)*(*this)((i+2)%3,(j+1)%3));
   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>
+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++) {
         t(i,j)+=a(i,k)*b(k,j);
       }
   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>
+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);
   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>
+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);
   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_>
+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>
+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>
+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>
+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>
+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) {
+T determinant(const TensorTyped<T,3,3>&t) {
   return t.determinant();
 }
 
+template <typename T>
 inline
-TensorGeneric<3,3> inverse(const TensorGeneric<3,3>&t) {
+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>
+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>
+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) {
+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) {
+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)) os<<t(i,j)<<" ";
       else os<<t(i,j);
@@ -485,6 +478,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
@@ -498,33 +494,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.
@@ -547,8 +550,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"
                               <<"Matrix:\n"<<mat<<"\n"
@@ -568,6 +571,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 //__PLUMED_tools_Tensor_h
+#endif
diff --git a/src/tools/Tools.h b/src/tools/Tools.h
index fee3059035..1b3cc78ebc 100644
--- a/src/tools/Tools.h
+++ b/src/tools/Tools.h
@@ -245,14 +245,14 @@ class Tools {
   static void set_to_zero(std::vector<VectorGeneric<n>> & vec) {
     unsigned s=vec.size();
     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();
diff --git a/src/tools/Vector.h b/src/tools/Vector.h
index 1a108ef25f..e79cff0751 100644
--- a/src/tools/Vector.h
+++ b/src/tools/Vector.h
@@ -74,6 +74,14 @@ int main(){
 
 */
 template<typename T, unsigned n> class VectorTyped;
+
+template<typename T, unsigned n>
+VectorTyped<T, n> delta(const VectorTyped<T, n>&,const VectorTyped<T, n>&);
+
+template<typename T, unsigned n>
+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>
@@ -134,11 +142,9 @@ class VectorTyped {
   template<typename U, typename J, unsigned m>
   friend VectorTyped<U, m> operator/(const VectorTyped<U, m>&,J);
 /// return v2-v1
-  template<typename U, unsigned m>
-  friend VectorTyped<U, m> delta(const VectorTyped<U, m>&v1,const VectorTyped<U, m>&v2);
+  friend VectorTyped delta<>(const VectorTyped&v1,const VectorTyped&v2);
 /// return v1 .scalar. v2
-  template<typename U, unsigned m>
-  friend U dotProduct(const VectorTyped<U, m>&,const VectorTyped<U, m>&);
+  friend 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
@@ -275,7 +281,7 @@ VectorTyped<T, n> operator*(VectorTyped<T, n> v,J s) {
 
 template<typename T, typename J, unsigned n>
 VectorTyped<T, n> operator/(const VectorTyped<T, n>&v,J s) {
-  return v*(1.0/s);
+  return v*(T(1.0)/s);
 }
 
 template<typename T, unsigned n>