diff --git a/adjointfield.hh b/adjointfield.hh
index 9580ef5..290eee4 100644
--- a/adjointfield.hh
+++ b/adjointfield.hh
@@ -222,7 +222,7 @@ 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) {
+template<typename Float, class Group>  adjointfield<Float, Group> operator*(const Float &x, const adjointfield<Float, Group> &A) {
   adjointfield<Float, Group> res(A.getLx(), A.getLy(), A.getLz(), A.getLt());
   for(size_t i = 0; i < A.getSize(); i++) {
     res[i] = x * A[i];
diff --git a/hmc-u1.cc b/hmc-u1.cc
index 1e8308a..b29f466 100644
--- a/hmc-u1.cc
+++ b/hmc-u1.cc
@@ -42,7 +42,7 @@ int main(int ac, char* av[]) {
     ("nsteps", po::value<size_t>(&n_steps)->default_value(1000), "n_steps")
     ("tau", po::value<double>(&tau)->default_value(1.), "trajectory length tau")
     ("exponent", po::value<size_t>(&exponent)->default_value(0), "exponent for rounding")
-    ("integrator", po::value<size_t>(&integs)->default_value(0), "itegration scheme to be used: 0=leapfrog, 1=lp_leapfrog, 2=omf4, 3=lp_omf4, 4=Euler, 5=RUTH, 6=omf2")
+    ("integrator", po::value<size_t>(&integs)->default_value(0), "itegration scheme to be used: 0=leapfrog, 1=lp_leapfrog, 2=omf4, 3=lp_omf4, 4=Euler, 5=RUTH, 6=omf2, 7=temperedLeapfrog")
     ;
 
   int err = parse_commandline(ac, av, desc, gparams);
@@ -71,7 +71,7 @@ int main(int ac, char* av[]) {
   cout << "## Plaquette after rnd trafo: " << plaquette*normalisation << endl; 
 
   // Molecular Dynamics parameters
-  md_params mdparams(n_steps, tau);
+  md_params<std::mt19937> mdparams(n_steps, tau);
   
   // generate list of monomials
   gaugemonomial<double, _u1> gm(0);
@@ -82,7 +82,7 @@ int main(int ac, char* av[]) {
   monomial_list.push_back(&gm);
   monomial_list.push_back(&km);
 
-  integrator<double, _u1> * md_integ = set_integrator<double, _u1>(integs, exponent);
+  integrator<double, _u1, std::mt19937> * md_integ = set_integrator<double, _u1, std::mt19937>(integs, exponent);
 
   std::ofstream os;
   if(gparams.icounter == 0) 
@@ -98,8 +98,9 @@ int main(int ac, char* av[]) {
     }
     // PRNG engine
     std::mt19937 engine(gparams.seed+i);
+    mdparams.setengine(&engine);
     // perform the MD update
-    md_update(U, engine, mdparams, monomial_list, *md_integ);
+    md_update(U, mdparams, monomial_list, *md_integ);
 
     double energy = gauge_energy(U);
     double E = 0., Q = 0.;
diff --git a/hmc.cc b/hmc.cc
index a6738eb..a82a26b 100644
--- a/hmc.cc
+++ b/hmc.cc
@@ -40,7 +40,7 @@ int main(int ac, char* av[]) {
     ("nsteps", po::value<size_t>(&n_steps)->default_value(1000), "n_steps")
     ("tau", po::value<double>(&tau)->default_value(1.), "trajectory length tau")
     ("exponent", po::value<size_t>(&exponent)->default_value(0), "exponent for rounding")
-    ("integrator", po::value<size_t>(&integs)->default_value(0), "itegration scheme to be used: 0=leapfrog, 1=lp_leapfrog, 2=omf4, 3=lp_omf4, 4=Euler, 5=RUTH, 6=omf2")
+    ("integrator", po::value<size_t>(&integs)->default_value(0), "itegration scheme to be used: 0=leapfrog, 1=lp_leapfrog, 2=omf4, 3=lp_omf4, 4=Euler, 5=RUTH, 6=omf2, 7=temperedLeapfrog")
     ;
 
   int err = parse_commandline(ac, av, desc, gparams);
@@ -69,7 +69,7 @@ int main(int ac, char* av[]) {
   cout << "## Plaquette after rnd trafo: " << plaquette*normalisation << endl; 
 
   // Molecular Dynamics parameters
-  md_params mdparams(n_steps, tau);
+  md_params<std::mt19937> mdparams(n_steps, tau);
   
   // generate list of monomials
   gaugemonomial<double, su2> gm(0);
@@ -80,7 +80,7 @@ int main(int ac, char* av[]) {
   monomial_list.push_back(&gm);
   monomial_list.push_back(&km);
 
-  integrator<double, su2> * md_integ = set_integrator<double, su2>(integs, exponent);
+  integrator<double, su2, std::mt19937> * md_integ = set_integrator<double, su2, std::mt19937>(integs, exponent);
 
   std::ofstream os;
   if(gparams.icounter == 0) 
@@ -96,8 +96,9 @@ int main(int ac, char* av[]) {
     }
     // PRNG engine
     std::mt19937 engine(gparams.seed+i);
+    mdparams.setengine(&engine);
     // perform the MD update
-    md_update(U, engine, mdparams, monomial_list, *md_integ);
+    md_update(U, mdparams, monomial_list, *md_integ);
 
     double energy = gauge_energy(U);
     rate += mdparams.getaccept();
diff --git a/integrator.hh b/integrator.hh
index a0218c8..ec29c67 100644
--- a/integrator.hh
+++ b/integrator.hh
@@ -10,23 +10,26 @@
 #include<list>
 #include<iostream>
 #include<cmath>
+#include<random>
 
-enum integrators { LEAPFROG = 0, LP_LEAPFROG = 1, OMF4 = 2, LP_OMF4 = 3, EULER = 4, RUTH = 5, OMF2 = 6};
+enum integrators { LEAPFROG = 0, LP_LEAPFROG = 1, OMF4 = 2, LP_OMF4 = 3, EULER = 4, RUTH = 5, OMF2 = 6, T_LEAPFROG=7};
 
 // virtual integrator class
-template<typename Float, class Group> class integrator{
+template<typename Float, class Group, class URNG> class integrator{
 public:
   integrator() {}
-  virtual void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
-                         md_params const &params, const bool restore = true) = 0;
+  virtual void integrate(std::list<monomial<Float, Group>*> &monomial_list,
+                         hamiltonian_field<Float, Group> &h, 
+                         md_params<URNG> &params, const bool restore = true) = 0;
 };
 
 // leapfrog integration scheme
-template<typename Float, class Group> class leapfrog : public integrator<Float, Group> {
+template<typename Float, class Group, class URNG> class leapfrog : public integrator<Float, Group, URNG> {
 public:
   leapfrog() {}
-  void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
-                 md_params const &params, const bool restore=true) {
+  void integrate(std::list<monomial<Float, Group>*> &monomial_list,
+                 hamiltonian_field<Float, Group> &h, 
+                 md_params<URNG> &params, const bool restore=true) {
 
     adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     
@@ -48,13 +51,57 @@ public:
   }
 };
 
+// tempered leapfrog integration scheme
+template<typename Float, class Group, class URNG> class tempered_leapfrog : public integrator<Float, Group, URNG> {
+public:
+  tempered_leapfrog() {}
+  void integrate(std::list<monomial<Float, Group>*> &monomial_list,
+                 hamiltonian_field<Float, Group> &h, 
+                 md_params<URNG> &params, const bool restore=true) {
+
+    const size_t nsteps = params.getnsteps();
+    //    std::normal_distribution<Float> normal(1., params.gettemperedwidth());
+    std::normal_distribution<Float> normal(1., params.gettemperedwidth());
+    Float a = Float(normal(*params.getengine()));
+    std::cout << "tempered a " << a << std::endl;
+    adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
+
+    bool is_even = (params.getnsteps()%2 == 0);
+    Float dtau = params.gettau()/Float(nsteps);
+    // initial half-step for the  momenta
+    deriv = sqrt(a) * deriv;
+    update_momenta(monomial_list, deriv, h, dtau/2.);
+    // first full step for gauge
+    update_gauge(h, dtau);
+    // nsteps-1 full steps
+    for(size_t i = 0; i < nsteps-1; i++) {
+      update_momenta(monomial_list, deriv, h, dtau/2);
+      if ((is_even && i < nsteps/2-1) || (!is_even && i < nsteps/2)) {
+        deriv = a * deriv;
+      }
+      if ((i >= nsteps/2)) {
+        deriv = (1./a) * deriv;
+      }
+      update_momenta(monomial_list, deriv, h, dtau/2);
+      update_gauge(h, dtau);
+    }
+
+    // final half-step for the momenta
+    update_momenta(monomial_list, deriv, h, dtau/2.);
+    deriv = (1./sqrt(a)) * deriv;
+    // restore SU
+    if(restore) h.U->restoreSU();
+  }
+};
+
+
 // leapfrog with rounding to allow for lower precision during MD update
-template<typename Float, class Group> class lp_leapfrog : public integrator<Float, Group> {
+template<typename Float, class Group, class URNG> class lp_leapfrog : public integrator<Float, Group, URNG> {
 public:
   lp_leapfrog() : n_prec(16) {}
   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) {
+                 md_params<URNG> &params, const bool restore = true) {
     adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     const size_t N = pow(10, n_prec);
 
@@ -79,11 +126,11 @@ private:
 };
 
 // OMF2 integration scheme (also called 2MN, second order minimal norm)
-template<typename Float, class Group> class omf2 : public integrator<Float, Group> {
+template<typename Float, class Group, class URNG> class omf2 : public integrator<Float, Group, URNG> {
 public:
   omf2() : lambda(0.1938) {}
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
-                 md_params const &params, const bool restore = true) {
+                 md_params<URNG> &params, const bool restore = true) {
 
     adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     Float dtau = params.gettau()/Float(params.getnsteps());
@@ -115,11 +162,11 @@ private:
 
 
 // OMF4 integration scheme
-template<typename Float, class Group> class omf4 : public integrator<Float, Group> {
+template<typename Float, class Group, class URNG> class omf4 : public integrator<Float, Group, URNG> {
 public:
   omf4() : rho(0.2539785108410595), theta(-0.03230286765269967), vartheta(0.08398315262876693), lambda(0.6822365335719091) {}
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
-                 md_params const &params, const bool restore = true) {
+                 md_params<URNG> &params, const bool restore = true) {
 
     adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     Float dtau = params.gettau()/Float(params.getnsteps());
@@ -154,12 +201,12 @@ private:
 };
 
 // OMF4 integration scheme in low precision
-template<typename Float, class Group> class lp_omf4 : public integrator<Float, Group> {
+template<typename Float, class Group, class URNG> class lp_omf4 : public integrator<Float, Group, URNG> {
 public:
   lp_omf4() : rho(0.2539785108410595), theta(-0.03230286765269967), vartheta(0.08398315262876693), lambda(0.6822365335719091), n_prec(16) {}
   lp_omf4(size_t n) : rho(0.2539785108410595), theta(-0.03230286765269967), vartheta(0.08398315262876693), lambda(0.6822365335719091), 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) {
+                 md_params<URNG> &params, const bool restore = true) {
 
     adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     const size_t N = pow(10, n_prec);
@@ -197,11 +244,11 @@ private:
 };
 
 // semi-implicit symplectic Euler integration scheme
-template<typename Float, class Group> class euler : public integrator<Float, Group> {
+template<typename Float, class Group, class URNG> class euler : public integrator<Float, Group, URNG> {
 public:
   euler() {}
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
-                 md_params const &params, const bool restore = true) {
+                 md_params<URNG> &params, const bool restore = true) {
 
     adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     Float dtau = params.gettau()/Float(params.getnsteps());
@@ -216,11 +263,11 @@ public:
 };
 
 // third order symplectic, but non-reversible integrator (Ruth)
-template<typename Float, class Group> class ruth : public integrator<Float, Group> {
+template<typename Float, class Group, class URNG> class ruth : public integrator<Float, Group, URNG> {
 public:
   ruth() {}
   void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h, 
-                 md_params const &params, const bool restore = true) {
+                 md_params<URNG> &params, const bool restore = true) {
 
     adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
     
@@ -240,39 +287,43 @@ public:
 };
 
 
-template<typename Float, class Group> integrator<Float, Group>* set_integrator(const size_t integs, const size_t exponent) {
-  integrator<Float, Group> * integ;
+template<typename Float, class Group, class URNG> integrator<Float, Group, URNG>* set_integrator(const size_t integs, const size_t exponent) {
+  integrator<Float, Group, URNG> * integ;
   if(static_cast<integrators>(integs) == LEAPFROG) {
-    integ = new leapfrog<Float, Group>();
+    integ = new leapfrog<Float, Group, URNG>();
     std::cerr << "leapfrog" << std::endl;
   }
   else if(static_cast<integrators>(integs) == LP_LEAPFROG) {
-    integ = new lp_leapfrog<Float, Group>(exponent);
+    integ = new lp_leapfrog<Float, Group, URNG>(exponent);
     std::cerr << "lp_leapfrog" << std::endl;
   }
   else if(static_cast<integrators>(integs) == OMF4) {
-    integ = new omf4<Float, Group>();
+    integ = new omf4<Float, Group, URNG>();
     std::cerr << "omf4" << std::endl;
   }
   else if(static_cast<integrators>(integs) == LP_OMF4) {
-    integ = new lp_omf4<Float, Group>(exponent);
+    integ = new lp_omf4<Float, Group, URNG>(exponent);
     std::cerr << "lp_omf4" << std::endl;
   }
   else if(static_cast<integrators>(integs) == EULER) {
-    integ = new euler<Float, Group>();
+    integ = new euler<Float, Group, URNG>();
     std::cerr << "euler" << std::endl;
   }
   else if(static_cast<integrators>(integs) == RUTH) {
-    integ = new ruth<Float, Group>();
+    integ = new ruth<Float, Group, URNG>();
     std::cerr << "ruth" << std::endl;
   }
   else if(static_cast<integrators>(integs) == OMF2) {
-    integ = new omf2<Float, Group>();
+    integ = new omf2<Float, Group, URNG>();
     std::cerr << "omf2" << std::endl;
   }
+  else if(static_cast<integrators>(integs) == T_LEAPFROG) {
+    integ = new tempered_leapfrog<Float, Group, URNG>();
+    std::cerr << "tempered leapfrog" << std::endl;
+  }
   else {
     std::cerr << "Integrator does not match, using default" << std::endl;
-    integ = new leapfrog<Float, Group>();
+    integ = new leapfrog<Float, Group, URNG>();
   }
   return integ;
 }
diff --git a/kramers.cc b/kramers.cc
index b9de2f0..f3f1a1f 100644
--- a/kramers.cc
+++ b/kramers.cc
@@ -61,7 +61,7 @@ int main(int ac, char* av[]) {
     hotstart(U, gparams.seed, gparams.heat);
   }
   // Molecular Dynamics parameters
-  md_params mdparams(n_steps, tau);
+  md_params<std::mt19937> mdparams(n_steps, tau);
   mdparams.setkmax(k_max);
   mdparams.setgamma(gamma);
 
@@ -81,7 +81,7 @@ int main(int ac, char* av[]) {
   monomial_list.push_back(&gm);
   monomial_list.push_back(&km);
 
-  integrator<double, su2> * md_integ = set_integrator<double, su2>(integs, exponent);
+  integrator<double, su2, std::mt19937> * md_integ = set_integrator<double, su2, std::mt19937>(integs, exponent);
 
   mdparams.setkmax(5);
   if(!gparams.acceptreject) mdparams.disableacceptreject();
@@ -98,8 +98,9 @@ int main(int ac, char* av[]) {
 
     // PRNG engine  
     std::mt19937 engine(gparams.seed + i);
+    mdparams.setengine(&engine);
     // perform the MD update
-    kramers_md_update(U, engine, mdparams, monomial_list, *md_integ);
+    kramers_md_update(U, mdparams, monomial_list, *md_integ);
 
     double energy = gauge_energy(U);
     rate += mdparams.getaccept();
diff --git a/kramers_md_update.hh b/kramers_md_update.hh
index d92b1cd..8236808 100644
--- a/kramers_md_update.hh
+++ b/kramers_md_update.hh
@@ -18,14 +18,14 @@ using std::vector;
 
 
 template<class URNG, typename Float, class Group> void kramers_md_update(gaugeconfig<su2> &U,
-                                                                         URNG &engine, 
-                                                                         md_params &params,
+                                                                         md_params<URNG> &params,
                                                                          std::list<monomial<Float, Group>*> &monomial_list, 
-                                                                         integrator<Float, Group> &md_integ) {
+                                                                         integrator<Float, Group, URNG> &md_integ) {
   adjointfield<Float, Group> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  URNG * engine = params.getengine();
   // generate standard normal distributed random momenta
   // normal distribution checked!
-  initnormal(engine, momenta);
+  initnormal(*engine, momenta);
 
   // for the accept reject step
   std::uniform_real_distribution<Float> uniform(0., 1.);
@@ -37,7 +37,7 @@ template<class URNG, typename Float, class Group> void kramers_md_update(gaugeco
 
   for(size_t k = 0; k < params.getkmax(); k++) {
     // first momenta update
-    initnormal<URNG, Float>(engine, eta);
+    initnormal<URNG, Float>(*engine, eta);
     Float epsilon = params.gettau()/Float(params.getnsteps());
     for(size_t i = 0; i < momenta.getSize(); i++) {
       momenta[i].seta(momenta[i].geta()*exp(-gamma*epsilon) + sqrt(1 - exp(-2*gamma*epsilon))*eta[i].geta());
@@ -77,7 +77,7 @@ template<class URNG, typename Float, class Group> void kramers_md_update(gaugeco
     params.setaccept(true);
     if(params.getacceptreject()) {
       if(delta_H > 0) {
-        if(uniform(engine) > exp(-delta_H)) {
+        if(uniform(*engine) > exp(-delta_H)) {
           params.setaccept(false);
         }
       }
diff --git a/lyapunov.hh b/lyapunov.hh
index dc7ff3b..030f769 100644
--- a/lyapunov.hh
+++ b/lyapunov.hh
@@ -22,9 +22,9 @@ using std::vector;
 
 template<typename Float, class URNG, class Group> void compute_lyapunov(gaugeconfig<Group> &U, 
                                                                         URNG &engine, 
-                                                                        md_params params,
+                                                                        md_params<URNG> params,
                                                                         std::list<monomial<Float, Group>*> &monomial_list, 
-                                                                        integrator<Float, Group> &md_integ, 
+                                                                        integrator<Float, Group, URNG> &md_integ, 
                                                                         std::string const &path, 
                                                                         const size_t d = 12) {
 
diff --git a/md_params.hh b/md_params.hh
index 086f6c8..df4f0da 100644
--- a/md_params.hh
+++ b/md_params.hh
@@ -1,12 +1,12 @@
 #pragma once
 
-class md_params {
+template<class URNG> class md_params {
 
 public:
   md_params(size_t _nsteps, double _tau, size_t n = 1000000000) : 
     nsteps(_nsteps), n_prec(n), kmax(5), exponent(16), tau(_tau), H(0.),
     deltaH(0.), deltadeltaH(0.), deltadeltaU(0.),
-    gamma(2.0),
+    gamma(2.0), omflambda(0.1938), temperedwidth(0.06),
     revtest(false), accept(true), acceptreject(true) {}
   size_t getnsteps() const {
     return nsteps;
@@ -93,6 +93,15 @@ public:
   size_t getexponent() const {
     return exponent;
   }
+  void setengine(URNG * _engine) {
+    engine = _engine;
+  }
+  URNG * getengine() {
+    return(engine);
+  }
+  double gettemperedwidth() {
+    return temperedwidth;
+  }
 private:
   size_t nsteps;
   size_t n_prec;
@@ -101,6 +110,9 @@ private:
   double tau;
   double H, deltaH, deltadeltaH, deltadeltaU;
   double gamma;
+  double omflambda;
+  double temperedwidth;
   bool revtest, accept, acceptreject;
+  URNG * engine;
 };
 
diff --git a/md_update.hh b/md_update.hh
index e4737d5..5d18003 100644
--- a/md_update.hh
+++ b/md_update.hh
@@ -16,14 +16,15 @@ using std::vector;
 
 
 template<class URNG, typename Float, class Group> void md_update(gaugeconfig<Group> &U,
-                                                                 URNG &engine, 
-                                                                 md_params &params,
+                                                                 md_params<URNG> &params,
                                                                  std::list<monomial<Float, Group>*> &monomial_list, 
-                                                                 integrator<Float, Group> &md_integ) {
+                                                                 integrator<Float, Group, URNG> &md_integ) {
   adjointfield<Float, Group> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
+  URNG * engine = params.getengine();
+
   // generate standard normal distributed random momenta
   // normal distribution checked!
-  initnormal(engine, momenta);
+  initnormal(*engine, momenta);
 
   std::uniform_real_distribution<Float> uniform(0., 1.);
 
@@ -53,7 +54,7 @@ template<class URNG, typename Float, class Group> void md_update(gaugeconfig<Gro
   // accept/reject step, if needed
   params.setaccept(true);
   if(delta_H > 0) {
-    if(uniform(engine) > exp(-delta_H)) {
+    if(uniform(*engine) > exp(-delta_H)) {
       params.setaccept(false);
     }
   }
diff --git a/measure-u1.cc b/measure-u1.cc
index 4f7f7b8..459358c 100644
--- a/measure-u1.cc
+++ b/measure-u1.cc
@@ -52,7 +52,6 @@ int main(int ac, char* av[]) {
     ("nsteps", po::value<size_t>(&n_steps)->default_value(1000), "n_steps")
     ("tau", po::value<double>(&tau)->default_value(1.), "trajectory length tau")
     ("exponent", po::value<size_t>(&exponent)->default_value(0), "exponent for rounding")
-    ("integrator", po::value<size_t>(&integs)->default_value(0), "itegration scheme to be used: 0=leapfrog, 1=lp_leapfrog, 2=omf4, 3=lp_omf4")
     ;
 
   int err = parse_commandline(ac, av, desc, gparams);
diff --git a/measure.cc b/measure.cc
index e83f278..f830d27 100644
--- a/measure.cc
+++ b/measure.cc
@@ -120,8 +120,8 @@ int main(int ac, char* av[]) {
       std::list<monomial<double, su2>*> monomial_list;
       monomial_list.push_back(&gm);
       monomial_list.push_back(&km);
-      md_params mdparams(n_steps, tau);
-      integrator<double, su2> * md_integ = set_integrator<double, su2>(integs, exponent);
+      md_params<std::mt19937> mdparams(n_steps, tau);
+      integrator<double, su2, std::mt19937> * md_integ = set_integrator<double, su2, std::mt19937>(integs, exponent);
 
       compute_lyapunov(U, engine, mdparams, monomial_list, *md_integ, os.str(), exponent);
       delete(md_integ);