Fix and compiles

Grid/qcd/hmc/GenericHMCrunner.h

@@ -140,17 +140,7 @@ private:
 
     // Can move this outside?
     typedef IntegratorType<SmearingPolicy> TheIntegrator;
-    // Metric
+    TheIntegrator MDynamics(UGrid, Parameters.MD, TheAction, Smearing);
-    //TrivialMetric<typename Implementation::Field> Mtr;
-    ConjugateGradient<LatticeGaugeField> CG(1.0e-8,10000);
-    LaplacianParams LapPar(0.0001, 1.0, 10000, 1e-8, 12, 64);
-//    RealD Kappa = 1.2;
-    RealD Kappa = Parameters.Kappa;
-    std::cout << GridLogMessage << "Kappa = " << Kappa << std::endl;
-
-    // Better to pass the generalised momenta to the integrator
-    LaplacianAdjointField<PeriodicGimplR> Laplacian(UGrid, CG, LapPar, Kappa);
-    TheIntegrator MDynamics(UGrid, Parameters.MD, TheAction, Smearing, Laplacian);
 
     if (Parameters.StartingType == "HotStart") {
       // Hot start

Grid/qcd/hmc/HMC.h

@@ -53,7 +53,6 @@ struct HMCparameters: Serializable {
                                   bool, MetropolisTest,
                                   Integer, NoMetropolisUntil,
                                   std::string, StartingType,
-				  RealD, Kappa,
                                   IntegratorParameters, MD)
 
   HMCparameters() {

Grid/qcd/hmc/integrators/Integrator.h

@@ -73,8 +73,7 @@ protected:
   double t_U;  // Track time passing on each level and for U and for P
   std::vector<double> t_P;  
 
-//  MomentaField P;
+  MomentaField P;
-  GeneralisedMomenta<FieldImplementation > P;
   SmearingPolicy& Smearer;
   RepresentationPolicy Representations;
   IntegratorParameters Params;
@@ -84,7 +83,7 @@ protected:
   void update_P(Field& U, int level, double ep) 
   {
     t_P[level] += ep;
-    update_P(P.Mom, U, level, ep);
+    update_P(P, U, level, ep);
 
     std::cout << GridLogIntegrator << "[" << level << "] P " << " dt " << ep << " : t_P " << t_P[level] << std::endl;
   }
@@ -112,21 +111,6 @@ protected:
     // input U actually not used in the fundamental case
     // Fundamental updates, include smearing
 
-    // Generalised momenta  
-    // Derivative of the kinetic term must be computed before
-    // Mom is the momenta and gets updated by the 
-    // actions derivatives
-    MomentaField MomDer(P.Mom.Grid());
-    P.M.ImportGauge(U);
-    P.DerivativeU(P.Mom, MomDer);
-    Mom -= MomDer * ep;
-
-    // Auxiliary fields
-    P.update_auxiliary_momenta(ep*0.5);
-    P.AuxiliaryFieldsDerivative(MomDer);
-    Mom -= MomDer * ep;
-    P.update_auxiliary_momenta(ep*0.5);
-
     for (int a = 0; a < as[level].actions.size(); ++a) {
       double start_full = usecond();
       Field force(U.Grid());
@@ -153,83 +137,9 @@ protected:
     as[level].apply(update_P_hireps, Representations, Mom, U, ep);
   }
 
-  void implicit_update_P(Field& U, int level, double ep, bool intermediate = false) {
-    t_P[level] += ep;
-
-    std::cout << GridLogIntegrator << "[" << level << "] P "
-              << " dt " << ep << " : t_P " << t_P[level] << std::endl;
-    // Fundamental updates, include smearing
-    MomentaField Msum(P.Mom.Grid());
-    Msum = Zero();
-    for (int a = 0; a < as[level].actions.size(); ++a) {
-      // Compute the force terms for the lagrangian part
-      // We need to compute the derivative of the actions
-      // only once
-      Field force(U.Grid());
-      conformable(U.Grid(), P.Mom.Grid());
-      Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
-      as[level].actions.at(a)->deriv(Us, force);  // deriv should NOT include Ta
-
-      std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl;
-      if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
-      force = FieldImplementation::projectForce(force);  // Ta for gauge fields
-      Real force_abs = std::sqrt(norm2(force) / U.Grid()->gSites());
-      std::cout << GridLogIntegrator << "|Force| site average: " << force_abs
-                << std::endl;
-      Msum += force;
-    }
-
-    MomentaField NewMom = P.Mom;
-    MomentaField OldMom = P.Mom;
-    double threshold = 1e-8;
-    P.M.ImportGauge(U);
-    MomentaField MomDer(P.Mom.Grid());
-    MomentaField MomDer1(P.Mom.Grid());
-    MomentaField AuxDer(P.Mom.Grid());
-    MomDer1 = Zero();
-    MomentaField diff(P.Mom.Grid());
-    double factor = 2.0;
-    if (intermediate){
-      P.DerivativeU(P.Mom, MomDer1);
-      factor = 1.0;
-    }
-
-    // Auxiliary fields
-    P.update_auxiliary_momenta(ep*0.5);
-    P.AuxiliaryFieldsDerivative(AuxDer);
-    Msum += AuxDer;
-    
-
-    // Here run recursively
-    int counter = 1;
-    RealD RelativeError;
-    do {
-      std::cout << GridLogIntegrator << "UpdateP implicit step "<< counter << std::endl;
-
-      // Compute the derivative of the kinetic term
-      // with respect to the gauge field
-      P.DerivativeU(NewMom, MomDer);
-      Real force_abs = std::sqrt(norm2(MomDer) / U.Grid()->gSites());
-      std::cout << GridLogIntegrator << "|Force| laplacian site average: " << force_abs
-                << std::endl;
-
-      NewMom = P.Mom - ep* 0.5 * (2.0*Msum + factor*MomDer + MomDer1);// simplify
-      diff = NewMom - OldMom;
-      counter++;
-      RelativeError = std::sqrt(norm2(diff))/std::sqrt(norm2(NewMom));
-      std::cout << GridLogIntegrator << "UpdateP RelativeError: " << RelativeError << std::endl;
-      OldMom = NewMom;
-    } while (RelativeError > threshold);
-
-    P.Mom = NewMom;
-
-    // update the auxiliary fields momenta    
-    P.update_auxiliary_momenta(ep*0.5);
-  }
-
   void update_U(Field& U, double ep) 
   {
-    update_U(P.Mom, U, ep);
+    update_U(P, U, ep);
 
     t_U += ep;
     int fl = levels - 1;
@@ -248,64 +158,15 @@ protected:
     Representations.update(U);  // void functions if fundamental representation
   }
 
-  void implicit_update_U(Field&U, double ep){
-    t_U += ep;
-    int fl = levels - 1;
-    std::cout << GridLogIntegrator << "   " << "[" << fl << "] U " << " dt " << ep << " : t_U " << t_U << std::endl;
-
-    MomentaField Mom1(P.Mom.Grid());
-    MomentaField Mom2(P.Mom.Grid());
-    RealD RelativeError;
-    Field diff(U.Grid());
-    Real threshold = 1e-8;
-    int counter = 1;
-    int MaxCounter = 100;
-
-    Field OldU = U;
-    Field NewU = U;
-
-    P.M.ImportGauge(U);
-    P.DerivativeP(Mom1); // first term in the derivative 
-
-    P.update_auxiliary_fields(ep*0.5);
-
-
-    MomentaField sum=Mom1;
-    do {
-      std::cout << GridLogIntegrator << "UpdateU implicit step "<< counter << std::endl;
-      
-      P.DerivativeP(Mom2); // second term in the derivative, on the updated U
-      sum = (Mom1 + Mom2);
-
-      for (int mu = 0; mu < Nd; mu++) {
-        auto Umu = PeekIndex<LorentzIndex>(U, mu);
-        auto Pmu = PeekIndex<LorentzIndex>(sum, mu);
-        Umu = expMat(Pmu, ep * 0.5, 12) * Umu;
-        PokeIndex<LorentzIndex>(NewU, ProjectOnGroup(Umu), mu);
-      }
-
-      diff = NewU - OldU;
-      RelativeError = std::sqrt(norm2(diff))/std::sqrt(norm2(NewU));
-      std::cout << GridLogIntegrator << "UpdateU RelativeError: " << RelativeError << std::endl;
-      
-      P.M.ImportGauge(NewU);
-      OldU = NewU; // some redundancy to be eliminated
-      counter++;
-    } while (RelativeError > threshold && counter < MaxCounter);
-
-    U = NewU;
-    P.update_auxiliary_fields(ep*0.5);
-  }
-
   virtual void step(Field& U, int level, int first, int last) = 0;
 
 public:
   Integrator(GridBase* grid, IntegratorParameters Par,
              ActionSet<Field, RepresentationPolicy>& Aset,
-             SmearingPolicy& Sm, Metric<MomentaField>& M)
+             SmearingPolicy& Sm)
     : Params(Par),
       as(Aset),
-      P(grid, M),
+      P(grid),
       levels(Aset.size()),
       Smearer(Sm),
       Representations(grid) 
@@ -342,9 +203,7 @@ public:
 
   void reverse_momenta()
   {
-//    P *= -1.0;
+    P *= -1.0;
-    P.Mom *= -1.0;
-    P.AuxMom *= -1.0;
   }
 
   // to be used by the actionlevel class to iterate
@@ -364,13 +223,10 @@ public:
   // Initialization of momenta and actions
   void refresh(Field& U, GridParallelRNG& pRNG) 
   {
-    assert(P.Mom.Grid() == U.Grid());
+    assert(P.Grid() == U.Grid());
     std::cout << GridLogIntegrator << "Integrator refresh\n";
 
-//    FieldImplementation::generate_momenta(P.Mom, pRNG);
+    FieldImplementation::generate_momenta(P, pRNG);
-    P.M.ImportGauge(U);
-    P.MomentaDistribution(pRNG);
-
 
     // Update the smeared fields, can be implemented as observer
     // necessary to keep the fields updated even after a reject
@@ -416,11 +272,9 @@ public:
 
     std::cout << GridLogIntegrator << "Integrator action\n";
 
-//    RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+    RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
-    P.M.ImportGauge(U);
+
-    RealD H = - P.MomentaAction();
     RealD Hterm;
-    std::cout << GridLogMessage << "Momentum action H_p = " << H << "\n";
 
     // Actions
     for (int level = 0; level < as.size(); ++level) {

Grid/qcd/hmc/integrators/Integrator_algorithm.h

@@ -101,8 +101,8 @@ public:
 
   std::string integrator_name(){return "LeapFrog";}
 
-  LeapFrog(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
+  LeapFrog(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
-    : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(grid, Par, Aset, Sm,M){};
+    : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(grid, Par, Aset, Sm){};
 
   void step(Field& U, int level, int _first, int _last) {
     int fl = this->as.size() - 1;
@@ -144,8 +144,8 @@ private:
 public:
   INHERIT_FIELD_TYPES(FieldImplementation);
 
-  MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
+  MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
-    : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(grid, Par, Aset, Sm,M){};
+    : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(grid, Par, Aset, Sm){};
 
   std::string integrator_name(){return "MininumNorm2";}
 
@@ -207,9 +207,9 @@ public:
   // Looks like dH scales as dt^4. tested wilson/wilson 2 level.
   ForceGradient(GridBase* grid, IntegratorParameters Par,
                 ActionSet<Field, RepresentationPolicy>& Aset,
-                SmearingPolicy& Sm, Metric<Field>& M)
+                SmearingPolicy& Sm)
     : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
-									    grid, Par, Aset, Sm,M){};
+									    grid, Par, Aset, Sm){};
 
   std::string integrator_name(){return "ForceGradient";}
   
@@ -271,139 +271,6 @@ public:
   }
 };
 
-////////////////////////////////
-// Riemannian Manifold HMC
-// Girolami et al
-////////////////////////////////
-
-
-
-// correct
-template <class FieldImplementation, class SmearingPolicy,
-          class RepresentationPolicy =
-              Representations<FundamentalRepresentation> >
-class ImplicitLeapFrog : public Integrator<FieldImplementation, SmearingPolicy,
-                                           RepresentationPolicy> {
- public:
-  typedef ImplicitLeapFrog<FieldImplementation, SmearingPolicy, RepresentationPolicy>
-      Algorithm;
-  INHERIT_FIELD_TYPES(FieldImplementation);
-
-  // Riemannian manifold metric operator
-  // Hermitian operator Fisher
-
-  std::string integrator_name(){return "ImplicitLeapFrog";}
-
-  ImplicitLeapFrog(GridBase* grid, IntegratorParameters Par,
-           ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
-      : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
-            grid, Par, Aset, Sm, M){};
-
-  void step(Field& U, int level, int _first, int _last) {
-    int fl = this->as.size() - 1;
-    // level  : current level
-    // fl     : final level
-    // eps    : current step size
-
-    // Get current level step size
-    RealD eps = this->Params.trajL/this->Params.MDsteps;
-    for (int l = 0; l <= level; ++l) eps /= this->as[l].multiplier;
-
-    int multiplier = this->as[level].multiplier;
-    for (int e = 0; e < multiplier; ++e) {
-      int first_step = _first && (e == 0);
-      int last_step = _last && (e == multiplier - 1);
-
-      if (first_step) {  // initial half step
-       this->implicit_update_P(U, level, eps / 2.0);
-      }
-
-      if (level == fl) {  // lowest level
-        this->implicit_update_U(U, eps);
-      } else {  // recursive function call
-        this->step(U, level + 1, first_step, last_step);
-      }
-
-      //int mm = last_step ? 1 : 2;
-      if (last_step){
-        this->update_P(U, level, eps / 2.0);
-      } else {
-      this->implicit_update_P(U, level, eps, true);// works intermediate step
-      //  this->update_P(U, level, eps); // looks not reversible
-      }
-    }
-  }
-};
-
-
-// This is not completely tested
-template <class FieldImplementation, class SmearingPolicy,
-          class RepresentationPolicy =
-              Representations<FundamentalRepresentation> >
-class ImplicitMinimumNorm2 : public Integrator<FieldImplementation, SmearingPolicy,
-                                       RepresentationPolicy> {
- private:
-  const RealD lambda = 0.1931833275037836;
-
- public:
-  INHERIT_FIELD_TYPES(FieldImplementation);
-
-  ImplicitMinimumNorm2(GridBase* grid, IntegratorParameters Par,
-               ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
-      : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
-            grid, Par, Aset, Sm, M){};
-
-  std::string integrator_name(){return "ImplicitMininumNorm2";}
-
-  void step(Field& U, int level, int _first, int _last) {
-    // level  : current level
-    // fl     : final level
-    // eps    : current step size
-
-    int fl = this->as.size() - 1;
-
-    RealD eps = this->Params.trajL/this->Params.MDsteps * 2.0;
-    for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier;
-
-    // Nesting:  2xupdate_U of size eps/2
-    // Next level is eps/2/multiplier
-
-    int multiplier = this->as[level].multiplier;
-    for (int e = 0; e < multiplier; ++e) {  // steps per step
-
-      int first_step = _first && (e == 0);
-      int last_step = _last && (e == multiplier - 1);
-
-      if (first_step) {  // initial half step
-        this->implicit_update_P(U, level, lambda * eps);
-      }
-
-      if (level == fl) {  // lowest level
-        this->implicit_update_U(U, 0.5 * eps);
-      } else {  // recursive function call
-        this->step(U, level + 1, first_step, 0);
-      }
-
-      this->implicit_update_P(U, level, (1.0 - 2.0 * lambda) * eps, true);
-
-      if (level == fl) {  // lowest level
-        this->implicit_update_U(U, 0.5 * eps);
-      } else {  // recursive function call
-        this->step(U, level + 1, 0, last_step);
-      }
-
-      //int mm = (last_step) ? 1 : 2;
-      //this->update_P(U, level, lambda * eps * mm);
-
-      if (last_step) {
-        this->update_P(U, level, eps * lambda);
-      } else {
-        this->implicit_update_P(U, level, lambda * eps*2.0, true);
-      }
-    }
-  }
-};
-
 NAMESPACE_END(Grid);
 
 #endif  // INTEGRATOR_INCLUDED

Grid/qcd/utils/A2Autils.h

@@ -1119,33 +1119,39 @@ void A2Autils<FImpl>::ContractFourQuarkColourDiagonal(const PropagatorField &WWV
   assert(gamma0.size()==gamma1.size());
   int Ng = gamma0.size();
 
+  // Make device accessible copy
+  Vector<Gamma> Gamma0v (Ng);
+  Vector<Gamma> Gamma1v (Ng);
+  Gamma *Gamma0 = & Gamma0v[0];
+  Gamma *Gamma1 = & Gamma1v[0];
+  for(int g=0;g<Ng;g++) {
+    Gamma0[g]=gamma0[g];
+    Gamma1[g]=gamma1[g];
+  }
+  
   GridBase *grid = WWVV0.Grid();
 
-  autoView(WWVV0_v , WWVV0,CpuRead);
-  autoView(WWVV1_v , WWVV1,CpuRead);
-  autoView(O_trtr_v, O_trtr,CpuWrite);
-  autoView(O_fig8_v, O_fig8,CpuWrite);
-  thread_for(ss,grid->oSites(),{
-
   typedef typename ComplexField::vector_object vobj;
+  autoView(WWVV0_v , WWVV0,AcceleratorRead);
+  autoView(WWVV1_v , WWVV1,AcceleratorRead);
+  autoView(O_trtr_v, O_trtr,AcceleratorWrite);
+  autoView(O_fig8_v, O_fig8,AcceleratorWrite);
+  accelerator_for(ss,grid->oSites(),vobj::Nsimd(),{
 
-    vobj v_trtr;
+    auto VV0 = WWVV0_v(ss);
-    vobj v_fig8;
+    auto VV1 = WWVV1_v(ss);
-
-    auto VV0 = WWVV0_v[ss];
-    auto VV1 = WWVV1_v[ss];
     
     for(int g=0;g<Ng;g++){
 
-      v_trtr = trace(VV0 * gamma0[g])* trace(VV1*gamma1[g]);
+      auto v_trtr = trace(VV0 * gamma0[g])* trace(VV1*gamma1[g]);
-      v_fig8 = trace(VV0 * gamma0[g] * VV1 * gamma1[g]);
+      auto v_fig8 = trace(VV0 * gamma0[g] * VV1 * gamma1[g]);
 
       if ( g==0 ) {
-	O_trtr_v[ss] = v_trtr; 
+	coalescedWrite(O_trtr_v[ss], v_trtr); 
-	O_fig8_v[ss] = v_fig8;
+	coalescedWrite(O_fig8_v[ss], v_fig8);
       } else { 
-	O_trtr_v[ss]+= v_trtr; 
+	coalescedWrite(O_trtr_v[ss], O_trtr_v(ss)+v_trtr); 
-	O_fig8_v[ss]+= v_fig8;
+	coalescedWrite(O_fig8_v[ss], O_fig8_v(ss)+v_fig8);
       }
       
     }
@@ -1165,25 +1171,36 @@ void A2Autils<FImpl>::ContractFourQuarkColourMix(const PropagatorField &WWVV0,
 
   GridBase *grid = WWVV0.Grid();
 
-  autoView( WWVV0_v , WWVV0,CpuRead);
+  // Make device accessible copy
-  autoView( WWVV1_v , WWVV1,CpuRead);
+  Vector<Gamma> Gamma0v (Ng);
-  autoView( O_trtr_v, O_trtr,CpuWrite);
+  Vector<Gamma> Gamma1v (Ng);
-  autoView( O_fig8_v, O_fig8,CpuWrite);
+  Gamma *Gamma0 = & Gamma0v[0];
+  Gamma *Gamma1 = & Gamma1v[0];
+  for(int g=0;g<Ng;g++) {
+    Gamma0[g]=gamma0[g];
+    Gamma1[g]=gamma1[g];
+  }
 
-  thread_for(ss,grid->oSites(),{
+  autoView( WWVV0_v , WWVV0,AcceleratorRead);
+  autoView( WWVV1_v , WWVV1,AcceleratorRead);
+  autoView( O_trtr_v, O_trtr,AcceleratorWrite);
+  autoView( O_fig8_v, O_fig8,AcceleratorWrite);
 
   typedef typename ComplexField::vector_object vobj;
+  accelerator_for(ss,grid->oSites(),vobj::Nsimd(),{
 
-    auto VV0 = WWVV0_v[ss];
+    auto VV0 = WWVV0_v(ss);
-    auto VV1 = WWVV1_v[ss];
+    auto VV1 = WWVV1_v(ss);
+
+    typedef decltype(trace(VV0)) scalar;
 
     for(int g=0;g<Ng;g++){
 
       auto VV0G = VV0 * gamma0[g];  // Spin multiply
       auto VV1G = VV1 * gamma1[g];
 
-      vobj v_trtr=Zero();
+      scalar v_trtr=Zero();
-      vobj v_fig8=Zero();
+      scalar v_fig8=Zero();
 
       /////////////////////////////////////////
       // Colour mixed
@@ -1234,11 +1251,11 @@ Bag [8,4]  fig8 (-227.58,3.58808e-17) trtr (-32.5776,1.83286e-17)     //  - 1602
       }}}}
 
       if ( g==0 ) {
-	O_trtr_v[ss] = v_trtr; 
+	coalescedWrite(O_trtr_v[ss] , v_trtr); 
-	O_fig8_v[ss] = v_fig8;
+	coalescedWrite(O_fig8_v[ss] , v_fig8);
       } else { 
-	O_trtr_v[ss]+= v_trtr; 
+	coalescedWrite(O_trtr_v[ss],O_trtr_v(ss) + v_trtr); 
-	O_fig8_v[ss]+= v_fig8;
+	coalescedWrite(O_fig8_v[ss],O_fig8_v(ss) + v_fig8);
       }
       
     }

tests/hmc/Test_hmc_WilsonGauge_Implicit.cc

@@ -1,155 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: neo <cossu@post.kek.jp>
-Author: Guido Cossu <guido.cossu@ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/Grid.h>
-
-namespace Grid{
-struct RMHMCActionParameters: Serializable {
-  GRID_SERIALIZABLE_CLASS_MEMBERS(RMHMCActionParameters,
-				  double, gauge_beta)
-
-  template <class ReaderClass >
-  RMHMCActionParameters(Reader<ReaderClass>& Reader){
-    read(Reader, "Action", *this);
-  }
-};
-}
-
-int main(int argc, char **argv) {
-  using namespace Grid;
-//  using namespace Grid::QCD;
-
-  Grid_init(&argc, &argv);
-  GridLogIntegrator.Active(1);
-  int threads = GridThread::GetThreads();
-  // here make a routine to print all the relevant information on the run
-  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
-
-   // Typedefs to simplify notation
-  typedef GenericHMCRunner<ImplicitMinimumNorm2> HMCWrapper;  // Uses the default minimum norm
-   // Serialiser
-//  typedef Grid::JSONReader       Serialiser;
-  typedef Grid::XmlReader       Serialiser;
-
-
-  HMCWrapper TheHMC;
-  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
-  
-  // Reader, file should come from command line
-  if (TheHMC.ParameterFile.empty()){
-    std::cout << "Input file not specified."
-              << "Use --ParameterFile option in the command line.\nAborting" 
-              << std::endl;
-    exit(1);
-  }
-  Serialiser Reader(TheHMC.ParameterFile);
-
-  RMHMCActionParameters ActionParams(Reader);  
-
-  // Grid from the command line
-  TheHMC.Resources.AddFourDimGrid("gauge");
- 
-
-  // Checkpointer definition
-  CheckpointerParameters CPparams(Reader);
-//  TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
-  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
-
-  RNGModuleParameters RNGpar(Reader);
-  TheHMC.Resources.SetRNGSeeds(RNGpar);
-
-  // Construct observables
-  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
-  typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
-  TheHMC.Resources.AddObservable<PlaqObs>();
-  TopologyObsParameters TopParams(Reader);
-  TheHMC.Resources.AddObservable<QObs>(TopParams);
-  /////////////////////////////////////////////////////////////
-  // Collect actions
-  WilsonGaugeActionR Waction(ActionParams.gauge_beta);
-
-  
-  ActionLevel<HMCWrapper::Field> Level1(1);
-  Level1.push_back(&Waction);
-  TheHMC.TheAction.push_back(Level1);
-  /////////////////////////////////////////////////////////////
-  TheHMC.Parameters.initialize(Reader);
-  TheHMC.Run(); 
-
-  Grid_finalize();
-
-} // main
-
-/* Examples for input files
-
-JSON
-
-{
-    "Checkpointer": {
-	"config_prefix": "ckpoint_json_lat",
-	"rng_prefix": "ckpoint_json_rng",
-	"saveInterval": 10,
-	"format": "IEEE64BIG"
-    },
-    "RandomNumberGenerator": {
-	"serial_seeds": "1 2 3 4 6",
-	"parallel_seeds": "55 7 8 9 11"
-    },
-    "Action":{
-	"gauge_beta": 5.8
-    },
-    "TopologyMeasurement":{
-	"interval": 1,
-	"do_smearing": true,
-	"Smearing":{
-	    "steps": 200,
-	    "step_size": 0.01,
-	    "meas_interval": 50,
-	    "maxTau": 2.0
-	}
-    },
-    "HMC":{
-	"StartTrajectory": 0,
-	"Trajectories": 10,
-	"MetropolisTest": true,
-	"NoMetropolisUntil": 10,
-	"StartingType": "HotStart",
-	"MD":{
-	    "name": "MinimumNorm2",
-	    "MDsteps": 40,
-	    "trajL": 1.0
-	}
-    }
-}
-
-
-XML example not provided yet
-
-*/