From 2d837fe543d8e4374e6025888a9d4d636b167c92 Mon Sep 17 00:00:00 2001
From: Carsten Urbach <curbach@gmx.de>
Date: Thu, 27 May 2021 23:08:03 +0200
Subject: [PATCH] refactoring of gaugeconfig and adjointfield, remove
 adjointfield next

---
 CMakeLists.txt       |  3 +-
 adjointfield.hh      | 21 +++++++-------
 exp_gauge.cc         | 11 --------
 exp_gauge.hh         | 13 ++++++++-
 gaugeconfig.hh       | 13 +++++++--
 gaugemonomial.hh     |  2 +-
 gradient_flow.hh     |  2 +-
 hamiltonian_field.hh |  4 +--
 integrator.hh        | 14 ++++-----
 kramers_md_update.hh |  6 ++--
 lyapunov.hh          |  4 +--
 md_update.hh         |  2 +-
 monomial.hh          |  4 +--
 site_types.hh        |  5 ++++
 su2.hh               | 67 ++++++++++++++++++++++++--------------------
 test.cc              |  3 ++
 update_momenta.hh    |  4 +--
 17 files changed, 99 insertions(+), 79 deletions(-)
 create mode 100644 site_types.hh

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 25d9b76..dc66794 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,7 +6,7 @@ cmake_minimum_required(VERSION 2.8.11)
 project(su2 C CXX)
 project(su2hmc C CXX)
 
-add_definitions(--std=c++11)
+add_definitions(--std=c++14)
 
 ## here we generate version.hh
 ## to track the git version of the executables
@@ -56,7 +56,6 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -D_USE_OMP_")
 # generated objects files or create a simple static library. The latter is done
 # here.
 add_library(su2
-    exp_gauge.cc
     print_program_options.cc
     parse_commandline.cc
 )
diff --git a/adjointfield.hh b/adjointfield.hh
index 9580ef5..5bb42cf 100644
--- a/adjointfield.hh
+++ b/adjointfield.hh
@@ -120,9 +120,10 @@ template<typename Float> struct adjoint_type<Float, _u1> {
 };
 
 
-template<typename Float, class Group> class adjointfield {
+template<typename Adjoint> class adjointfield {
 public:
-  using value_type = typename adjoint_type<Float, Group>::type;
+  //using value_type = typename adjoint_type<Float, Group>::type;
+  using value_type = Adjoint;
   adjointfield(const size_t Lx, const size_t Ly, const size_t Lz, const size_t Lt, const size_t ndims = 4) : 
     Lx(Lx), Ly(Ly), Lz(Lz), Lt(Lt), volume(Lx*Ly*Lz*Lt), ndims(ndims) {
     data.resize(volume*4);
@@ -161,7 +162,7 @@ public:
   size_t getSize() const {
     return(volume*ndims);
   }
-  void operator=(const adjointfield<Float, Group> &U) {
+  void operator=(const adjointfield<Adjoint> &U) {
     Lx = U.getLx();
     Ly = U.getLz();
     Lz = U.getLz();
@@ -222,15 +223,15 @@ template<typename Float> inline adjointu1<Float> operator*(const Float &x, const
 }
 
 
-template<typename Float, class Group>  adjointfield<Float, su2> operator*(const Float &x, const adjointfield<Float, Group> &A) {
-  adjointfield<Float, Group> res(A.getLx(), A.getLy(), A.getLz(), A.getLt());
+template<typename Float, typename Adjoint>  adjointfield<Adjoint> operator*(const Float &x, const adjointfield<Adjoint> &A) {
+  adjointfield<Adjoint> res(A.getLx(), A.getLy(), A.getLz(), A.getLt());
   for(size_t i = 0; i < A.getSize(); i++) {
     res[i] = x * A[i];
   }
   return res;
 }
 
-template<typename Float> Float operator*(const adjointfield<Float, su2> &A, const adjointfield<Float, su2> &B) {
+template<typename Float> Float operator*(const adjointfield<adjointsu2<Float>> &A, const adjointfield<adjointsu2<Float>> &B) {
   Float res = 0.;
   assert(A.getSize() == B.getSize());
   for(size_t i = 0; i < A.getSize(); i++) {
@@ -239,7 +240,7 @@ template<typename Float> Float operator*(const adjointfield<Float, su2> &A, cons
   return res;
 }
 
-template<typename Float> Float operator*(const adjointfield<Float, _u1> &A, const adjointfield<Float, _u1> &B) {
+template<typename Float> Float operator*(const adjointfield<adjointu1<Float>> &A, const adjointfield<adjointu1<Float>> &B) {
   Float res = 0.;
   assert(A.getSize() == B.getSize());
   for(size_t i = 0; i < A.getSize(); i++) {
@@ -249,7 +250,7 @@ template<typename Float> Float operator*(const adjointfield<Float, _u1> &A, cons
 }
 
 
-template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<Float, su2> &A) {
+template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<adjointsu2<Float>> &A) {
   std::normal_distribution<double> normal(0., 1.);
   for(size_t i = 0; i < A.getSize(); i++) {
     A[i].seta(Float(normal(engine)));
@@ -259,7 +260,7 @@ template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<
   return;
 }
 
-template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<Float, _u1> &A) {
+template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<adjointu1<Float>> &A) {
   std::normal_distribution<double> normal(0., 1.);
   for(size_t i = 0; i < A.getSize(); i++) {
     A[i].seta(Float(normal(engine)));
@@ -267,7 +268,7 @@ template<class URNG, typename Float> void initnormal(URNG &engine, adjointfield<
   return;
 }
 
-template<typename Float, class Group> inline void zeroadjointfield(adjointfield<Float, Group> &A) {
+template<typename Adjoint> inline void zeroadjointfield(adjointfield<Adjoint> &A) {
   for(size_t i = 0; i < A.getSize(); i++) {
     A[i].setzero();
   }
diff --git a/exp_gauge.cc b/exp_gauge.cc
index 3708758..e03d58c 100644
--- a/exp_gauge.cc
+++ b/exp_gauge.cc
@@ -5,16 +5,5 @@
 #include<vector>
 #include<complex>
 
-_su2 exp(adjointsu2<double> const & x) {
-  double a = x.geta(), b = x.getb(), c = x.getc();
-  // normalise the vector
-  const double alpha = sqrt(a*a+b*b+c*c);
-  std::vector<double> n = {a/alpha, b/alpha, c/alpha};
-
-  const double salpha = sin(alpha);
-  _su2 res(Complex(cos(alpha), salpha*n[2]), 
-           salpha*Complex(n[1], n[0]));
-  return res;
-}
 
 
diff --git a/exp_gauge.hh b/exp_gauge.hh
index 68231f9..e1b2462 100644
--- a/exp_gauge.hh
+++ b/exp_gauge.hh
@@ -4,7 +4,18 @@
 #include "u1.hh"
 #include "adjointfield.hh"
 
-_su2 exp(adjointsu2<double> const & x);
+template<typename Float=double>
+_su2<Float> exp(adjointsu2<Float> const & x) {
+  double a = x.geta(), b = x.getb(), c = x.getc();
+  // normalise the vector
+  const double alpha = sqrt(a*a+b*b+c*c);
+  std::vector<double> n = {a/alpha, b/alpha, c/alpha};
+
+  const double salpha = sin(alpha);
+  _su2<Float> res(typename _su2<Float>::Complex(cos(alpha), salpha*n[2]), 
+                  salpha*(typename _su2<Float>::Complex(n[1], n[0])));
+  return res;
+}
 
 inline _u1 exp(adjointu1<double> const & x) {
   return _u1(x.geta());
diff --git a/gaugeconfig.hh b/gaugeconfig.hh
index c1f1f8a..20d1b0e 100644
--- a/gaugeconfig.hh
+++ b/gaugeconfig.hh
@@ -1,7 +1,6 @@
 #pragma once
 
-#include"su2.hh"
-#include"u1.hh"
+#include"site_types.hh"
 #include"random_element.hh"
 #include<random>
 #include<vector>
@@ -10,7 +9,7 @@
 #include<complex>
 #include<iostream>
 #include<cassert>
-
+#include<typeinfo>
 using std::vector;
 
 template<class T> class gaugeconfig {
@@ -68,6 +67,14 @@ public:
     }
   }
 
+  //  template<class Q=T>
+  template<class Q = T, typename = std::enable_if_t< (typeid(Q) == typeid(su2)) > >
+  void flipsign() {
+    for(size_t i = 0; i < getSize(); i++) {
+      data[i].flipsign();
+    }
+  }
+ 
   void operator=(const gaugeconfig &U) {
     volume = U.getVolume();
     data.resize(U.getSize());
diff --git a/gaugemonomial.hh b/gaugemonomial.hh
index 8703c75..ccab459 100644
--- a/gaugemonomial.hh
+++ b/gaugemonomial.hh
@@ -24,7 +24,7 @@ public:
     monomial<Float, Group>::Hnew = h.U->getBeta()*(h.U->getVolume()*6 - gauge_energy(*(h.U))/double(h.U->getNc()));
     return;
   }
-  void derivative(adjointfield<Float, Group> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac = 1.) const override {
+  void derivative(adjointfield<typename adjoint_type<Float, Group>::type> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac = 1.) const override {
     std::vector<size_t> x = {0, 0, 0, 0};
     typedef typename accum_type<Group>::type accum;
 #pragma omp parallel for
diff --git a/gradient_flow.hh b/gradient_flow.hh
index 0a3a320..5e62e17 100644
--- a/gradient_flow.hh
+++ b/gradient_flow.hh
@@ -63,7 +63,7 @@ template<class Group> void gradient_flow(gaugeconfig<Group> &U, std::string cons
   E[2] = density;
 
   gaugeconfig<Group> Vt(U);
-  adjointfield<double, Group> deriv(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<double, Group>::type> deriv(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   hamiltonian_field<double, Group> h(deriv, Vt);
   gaugemonomial<double, Group> SW(0);
 
diff --git a/hamiltonian_field.hh b/hamiltonian_field.hh
index aa2546f..e7c2530 100644
--- a/hamiltonian_field.hh
+++ b/hamiltonian_field.hh
@@ -5,8 +5,8 @@
 #include<vector>
 
 template<typename Float, class Group> struct hamiltonian_field {
-  adjointfield<Float, Group> * momenta;
+  adjointfield<typename adjoint_type<Float, Group>::type> * momenta;
   gaugeconfig<Group> * U;
-  hamiltonian_field(adjointfield<Float, Group> &momenta, gaugeconfig<Group> &U) :
+  hamiltonian_field(adjointfield<typename adjoint_type<Float, Group>::type> &momenta, gaugeconfig<Group> &U) :
     momenta(&momenta), U(&U) {}
 };
diff --git a/integrator.hh b/integrator.hh
index a0218c8..c7cdf45 100644
--- a/integrator.hh
+++ b/integrator.hh
@@ -28,7 +28,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore=true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     
     Float dtau = params.gettau()/Float(params.getnsteps());
     // initial half-step for the  momenta
@@ -55,7 +55,7 @@ public:
   lp_leapfrog(size_t n) : n_prec(n) {}
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     const size_t N = pow(10, n_prec);
 
     Float dtau = params.gettau()/Float(params.getnsteps());
@@ -85,7 +85,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     Float dtau = params.gettau()/Float(params.getnsteps());
     Float oneminus2lambda = (1.-2.*lambda);
 
@@ -121,7 +121,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     Float dtau = params.gettau()/Float(params.getnsteps());
     Float eps[10] = {rho*dtau, lambda*dtau, 
                  theta*dtau, 0.5*(1-2.*(lambda+vartheta))*dtau, 
@@ -161,7 +161,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     const size_t N = pow(10, n_prec);
 
     Float dtau = params.gettau()/Float(params.getnsteps());
@@ -203,7 +203,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     Float dtau = params.gettau()/Float(params.getnsteps());
     // nsteps full steps
     for(size_t i = 0; i < params.getnsteps(); i++) {
@@ -222,7 +222,7 @@ public:
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
                  md_params const &params, const bool restore = true) {
 
-    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+    adjointfield<typename adjoint_type<Float, Group>::type> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     
     Float dtau = params.gettau()/Float(params.getnsteps());
     // nsteps full steps
diff --git a/kramers_md_update.hh b/kramers_md_update.hh
index d92b1cd..5dec331 100644
--- a/kramers_md_update.hh
+++ b/kramers_md_update.hh
@@ -22,7 +22,7 @@ template<class URNG, typename Float, class Group> void kramers_md_update(gaugeco
                                                                          md_params &params,
                                                                          std::list<monomial<Float, Group>*> &monomial_list, 
                                                                          integrator<Float, Group> &md_integ) {
-  adjointfield<Float, Group> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   // generate standard normal distributed random momenta
   // normal distribution checked!
   initnormal(engine, momenta);
@@ -31,8 +31,8 @@ template<class URNG, typename Float, class Group> void kramers_md_update(gaugeco
   std::uniform_real_distribution<Float> uniform(0., 1.);
 
   const double gamma = params.getgamma();
-  adjointfield<Float, Group> eta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
-  adjointfield<Float, Group> momenta_old(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> eta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta_old(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   gaugeconfig<su2> U_old(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims(), U.getBeta());
 
   for(size_t k = 0; k < params.getkmax(); k++) {
diff --git a/lyapunov.hh b/lyapunov.hh
index dc7ff3b..36e021f 100644
--- a/lyapunov.hh
+++ b/lyapunov.hh
@@ -30,10 +30,10 @@ template<typename Float, class URNG, class Group> void compute_lyapunov(gaugecon
 
   const size_t n = pow(10, d);
 
-  adjointfield<Float, Group> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   // generate standard normal distributed random momenta
   initnormal(engine, momenta);
-  adjointfield<Float, Group> momenta2(momenta);
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta2(momenta);
   
   // generate copy of U, but round to d decimal digits
   gaugeconfig<Group> U2(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getBeta());
diff --git a/md_update.hh b/md_update.hh
index e4737d5..d3d4207 100644
--- a/md_update.hh
+++ b/md_update.hh
@@ -20,7 +20,7 @@ template<class URNG, typename Float, class Group> void md_update(gaugeconfig<Gro
                                                                  md_params &params,
                                                                  std::list<monomial<Float, Group>*> &monomial_list, 
                                                                  integrator<Float, Group> &md_integ) {
-  adjointfield<Float, Group> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  adjointfield<typename adjoint_type<Float, Group>::type> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
   // generate standard normal distributed random momenta
   // normal distribution checked!
   initnormal(engine, momenta);
diff --git a/monomial.hh b/monomial.hh
index 97df0e2..c8ae1e6 100644
--- a/monomial.hh
+++ b/monomial.hh
@@ -10,7 +10,7 @@ public:
   monomial() : Hold(0.), Hnew(0.), timescale(0), mdactive(true) {}
   virtual void heatbath(hamiltonian_field<Float, Group> const &h) = 0;
   virtual void accept(hamiltonian_field<Float, Group> const &h) = 0;
-  virtual void derivative(adjointfield<Float, Group> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac) const = 0;
+  virtual void derivative(adjointfield<typename adjoint_type<Float, Group>::type> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac) const = 0;
   void reset() {
     Hold = 0.;
     Hnew = 0.;
@@ -47,7 +47,7 @@ public:
   void accept(hamiltonian_field<Float, Group> const &h) override {
     monomial<Float, Group>::Hnew = 0.5 * ((*h.momenta)*(*h.momenta));
   }
-  virtual void derivative(adjointfield<Float, Group> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac) const {}
+  virtual void derivative(adjointfield<typename adjoint_type<Float, Group>::type> &deriv, hamiltonian_field<Float, Group> const &h, const Float fac) const {}
 };
 
 #include"gaugemonomial.hh"
diff --git a/site_types.hh b/site_types.hh
new file mode 100644
index 0000000..07a2a92
--- /dev/null
+++ b/site_types.hh
@@ -0,0 +1,5 @@
+#pragma once
+
+#include"su2.hh"
+#include"u1.hh"
+#include"adjointfield.hh"
diff --git a/su2.hh b/su2.hh
index 75dbaeb..1ca2d6b 100644
--- a/su2.hh
+++ b/su2.hh
@@ -3,31 +3,32 @@
 #include"accum_type.hh"
 #include"traceless_antiherm.hh"
 #include<complex>
-//#include<cmath>
 
-using Complex = std::complex<double>;
 
-
-class _su2 {
+template<typename Float=double> class _su2 {
 public:
+  typedef std::complex<Float> Complex;
   const size_t N_c = 2;
   explicit _su2() : a(0), b(0) {}
   explicit _su2(Complex a, Complex b) : a(a), b(b) {}
-  _su2(const _su2& U) : a(U.a), b(U.b) {}
-
-  friend inline _su2 operator+(const _su2 &U1, const _su2 &U2);
-  friend inline _su2 operator-(const _su2 &U1, const _su2 &U2);
-  friend inline _su2 operator*(const _su2 &U1, const _su2 &U2);
-  _su2& operator*=(const _su2 &U1) {
+  _su2(const _su2<Float>& U) : a(U.a), b(U.b) {}
+
+  template<typename F>
+  friend inline _su2<F> operator+(const _su2<F> &U1, const _su2<F> &U2);
+  template<typename F>
+  friend inline _su2<F> operator-(const _su2<F> &U1, const _su2<F> &U2);
+  template<typename F>
+  friend inline _su2<F> operator*(const _su2<F> &U1, const _su2<F> &U2);
+  _su2<Float>& operator*=(const _su2<Float> &U1) {
     Complex a = this->a;
     this->a = a*U1.a - this->b*std::conj(U1.b);
     this->b = a*U1.b + this->b*std::conj(U1.a);
     return *this;
   }
-  _su2 round(size_t n) const {
+  _su2<Float> round(size_t n) const {
     double dn = n;
-    return _su2(Complex(std::round(std::real(a)*dn), std::round(std::imag(a)*dn))/dn,
-                Complex(std::round(std::real(b)*dn), std::round(std::imag(b)*dn))/dn);
+    return _su2<Float>(Complex(std::round(std::real(a)*dn), std::round(std::imag(a)*dn))/dn,
+                       Complex(std::round(std::real(b)*dn), std::round(std::imag(b)*dn))/dn);
   }
 
   inline Complex geta() const {
@@ -36,11 +37,11 @@ public:
   inline Complex getb() const {
     return(b);
   }
-  inline void operator=(const _su2 &U) {
+  inline void operator=(const _su2<Float> &U) {
     a = U.geta();
     b = U.getb();
   }
-  inline void operator+=(const _su2 &U) {
+  inline void operator+=(const _su2<Float> &U) {
     a += U.a;
     b += U.b;
   }
@@ -48,8 +49,8 @@ public:
     a = _a;
     b = _b;
   }
-  inline _su2 dagger() const {
-    return(_su2(std::conj(a), -b));
+  inline _su2<Float> dagger() const {
+    return(_su2<Float>(std::conj(a), -b));
   }
   inline double retrace() {
     return(2.*std::real(a));
@@ -67,42 +68,46 @@ private:
   Complex a, b;
 };
 
-inline double retrace(_su2 const &U) {
+template<typename Float=double>
+inline double retrace(_su2<Float> const &U) {
   double a = std::real(U.geta());
   return(2*a);
 }
-
-inline Complex trace(_su2 const &U) {
+template<typename Float=double>
+inline typename _su2<Float>::Complex trace(_su2<Float> const &U) {
   double a = std::real(U.geta());
-  return(Complex(2*a, 0.));
+  return(typename _su2<Float>::Complex(2*a, 0.));
 }
 
-
-template<> inline _su2 traceless_antiherm(const _su2& x) {
-  return(_su2(0.5*(x.geta()-std::conj(x.geta())), x.getb()));
+template<typename Float=double>
+inline _su2<Float> traceless_antiherm(const _su2<Float>& x) {
+  return(_su2<Float>(0.5*(x.geta()-std::conj(x.geta())), x.getb()));
 }
 
-inline _su2 operator*(const _su2 &U1, const _su2 &U2) {
-  _su2 res;
+template<typename Float=double>
+inline _su2<Float> operator*(const _su2<Float> &U1, const _su2<Float> &U2) {
+  _su2<Float> res;
   res.a = U1.a*U2.a - U1.b*std::conj(U2.b);
   res.b = U1.a*U2.b + U1.b*std::conj(U2.a);
   return(res);
 }
 
-inline _su2 operator+(const _su2 &U1, const _su2 &U2) {
-  _su2 res;
+template<typename Float=double>
+inline _su2<Float> operator+(const _su2<Float> &U1, const _su2<Float> &U2) {
+  _su2<Float> res;
   res.a = U1.a + U2.a;
   res.b = U1.b + U2.b;
   return(res);
 }
 
-inline _su2 operator-(const _su2 &U1, const _su2 &U2) {
-  _su2 res;
+template<typename Float=double>
+inline _su2<Float> operator-(const _su2<Float> &U1, const _su2<Float> &U2) {
+  _su2<Float> res;
   res.a = U1.a - U2.a;
   res.b = U1.b - U2.b;
   return(res);
 }
 
 
-using su2 = _su2;
+using su2 = _su2<double>;
 
diff --git a/test.cc b/test.cc
index 8f6cb47..4bf164a 100644
--- a/test.cc
+++ b/test.cc
@@ -9,6 +9,7 @@
 
 #include<iostream>
 #include<vector>
+#include<typeinfo>
 
 using std::vector;
 using std::cout;
@@ -112,5 +113,7 @@ int main() {
   Q = 0;
   energy_density(cU, res, Q);
   cout << "Charge after random gauge trafo: " << Q << endl;
+
+  cout << typeid(U).name() << endl;
   return(0);
 }
diff --git a/update_momenta.hh b/update_momenta.hh
index 91ca93d..9ca74c8 100644
--- a/update_momenta.hh
+++ b/update_momenta.hh
@@ -8,7 +8,7 @@
 #include<iostream>
 
 template<typename Float, class Group> void update_momenta(std::list<monomial<Float, Group>*> &monomial_list, 
-                                                          adjointfield<Float, Group> &deriv, hamiltonian_field<Float, Group> &h, 
+                                                          adjointfield<typename adjoint_type<Float, Group>::type> &deriv, hamiltonian_field<Float, Group> &h, 
                                                           const double dtau) {
 
   zeroadjointfield(deriv);
@@ -33,7 +33,7 @@ template<typename Float, class Group> void update_momenta(std::list<monomial<Flo
 }
 
 template<typename Float, class Group> void round_and_update_momenta(std::list<monomial<Float, Group>*> &monomial_list, 
-                                                                    adjointfield<Float, Group> &deriv, hamiltonian_field<Float, Group> &h, 
+                                                                    adjointfield<typename adjoint_type<Float, Group>::type> &deriv, hamiltonian_field<Float, Group> &h, 
                                                                     const double dtau, const size_t n) {
   
   zeroadjointfield(deriv);