Fix for the CI

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -164,6 +164,7 @@ namespace QCD {
     public:
 
     static const int Dimension = Representation::Dimension;
+    static const bool isFundamental = Representation::isFundamental;
     static const bool LsVectorised=false;
     static const int Nhcs = Options::Nhcs;
 
@@ -298,27 +299,28 @@ namespace QCD {
   ////////////////////////////////////////////////////////////////////////////////////
   // Single flavour four spinors with colour index, 5d redblack
   ////////////////////////////////////////////////////////////////////////////////////
-template<class S,int Nrepresentation=Nc, class Options=CoeffReal>
+template<class S,class Representation = FundamentalRepresentation, class Options=CoeffReal>
-class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > { 
+class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Representation::Dimension> > { 
   public:
 
-  typedef PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;
+  typedef PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension> > Gimpl;
   INHERIT_GIMPL_TYPES(Gimpl);
 
-  static const int Dimension = Nrepresentation;
+  static const int Dimension = Representation::Dimension;
+  static const bool isFundamental = Representation::isFundamental;
   static const bool LsVectorised=true;
   static const int Nhcs = Options::Nhcs;
       
   typedef typename Options::_Coeff_t Coeff_t;      
   typedef typename Options::template PrecisionMapper<Simd>::LowerPrecVector SimdL;
   
-  template <typename vtype> using iImplSpinor            = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
+  template <typename vtype> using iImplSpinor            = iScalar<iVector<iVector<vtype, Dimension>, Ns> >;
-  template <typename vtype> using iImplPropagator        = iScalar<iMatrix<iMatrix<vtype, Nrepresentation>, Ns> >;
+  template <typename vtype> using iImplPropagator        = iScalar<iMatrix<iMatrix<vtype, Dimension>, Ns> >;
-  template <typename vtype> using iImplHalfSpinor        = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >;
+  template <typename vtype> using iImplHalfSpinor        = iScalar<iVector<iVector<vtype, Dimension>, Nhs> >;
-  template <typename vtype> using iImplHalfCommSpinor    = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhcs> >;
+  template <typename vtype> using iImplHalfCommSpinor    = iScalar<iVector<iVector<vtype, Dimension>, Nhcs> >;
-  template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>;
+  template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
-  template <typename vtype> using iImplGaugeField        = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd>;
+  template <typename vtype> using iImplGaugeField        = iVector<iScalar<iMatrix<vtype, Dimension> >, Nd>;
-  template <typename vtype> using iImplGaugeLink         = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
+  template <typename vtype> using iImplGaugeLink         = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
   
   typedef iImplSpinor<Simd>            SiteSpinor;
   typedef iImplPropagator<Simd>        SitePropagator;
@@ -354,8 +356,8 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
                        const SiteHalfSpinor &chi, int mu, StencilEntry *SE,
                        StencilImpl &St) {
     SiteGaugeLink UU;
-    for (int i = 0; i < Nrepresentation; i++) {
+    for (int i = 0; i < Dimension; i++) {
-      for (int j = 0; j < Nrepresentation; j++) {
+      for (int j = 0; j < Dimension; j++) {
         vsplat(UU()()(i, j), U(mu)()(i, j));
       }
     }
@@ -367,8 +369,8 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
                            const SitePropagator &chi,
                            int mu) {
     SiteGaugeLink UU;
-    for (int i = 0; i < Nrepresentation; i++) {
+    for (int i = 0; i < Dimension; i++) {
-      for (int j = 0; j < Nrepresentation; j++) {
+      for (int j = 0; j < Dimension; j++) {
         vsplat(UU()()(i, j), U(mu)()(i, j));
       }
     }
@@ -472,25 +474,26 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
     ////////////////////////////////////////////////////////////////////////////////////////
     // Flavour doubled spinors; is Gparity the only? what about C*?
     ////////////////////////////////////////////////////////////////////////////////////////
-template <class S, int Nrepresentation, class Options=CoeffReal>
+template <class S, class Representation = FundamentalRepresentation, class Options=CoeffReal>
-class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
+class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Representation::Dimension> > {
  public:
 
- static const int Dimension = Nrepresentation;
+ static const int Dimension = Representation::Dimension;
+ static const bool isFundamental = Representation::isFundamental;
  static const int Nhcs = Options::Nhcs;
  static const bool LsVectorised=false;
 
- typedef ConjugateGaugeImpl< GaugeImplTypes<S,Nrepresentation> > Gimpl;
+ typedef ConjugateGaugeImpl< GaugeImplTypes<S,Dimension> > Gimpl;
  INHERIT_GIMPL_TYPES(Gimpl);
 
  typedef typename Options::_Coeff_t Coeff_t;
  typedef typename Options::template PrecisionMapper<Simd>::LowerPrecVector SimdL;
       
- template <typename vtype> using iImplSpinor            = iVector<iVector<iVector<vtype, Nrepresentation>, Ns>,   Ngp>;
+ template <typename vtype> using iImplSpinor            = iVector<iVector<iVector<vtype, Dimension>, Ns>,   Ngp>;
- template <typename vtype> using iImplPropagator        = iVector<iMatrix<iMatrix<vtype, Nrepresentation>, Ns>,   Ngp>;
+ template <typename vtype> using iImplPropagator        = iVector<iMatrix<iMatrix<vtype, Dimension>, Ns>,   Ngp>;
- template <typename vtype> using iImplHalfSpinor        = iVector<iVector<iVector<vtype, Nrepresentation>, Nhs>,  Ngp>;
+ template <typename vtype> using iImplHalfSpinor        = iVector<iVector<iVector<vtype, Dimension>, Nhs>,  Ngp>;
- template <typename vtype> using iImplHalfCommSpinor    = iVector<iVector<iVector<vtype, Nrepresentation>, Nhcs>, Ngp>;
+ template <typename vtype> using iImplHalfCommSpinor    = iVector<iVector<iVector<vtype, Dimension>, Nhcs>, Ngp>;
- template <typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>, Ngp>;
+ template <typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>, Ngp>;
 
  typedef iImplSpinor<Simd>            SiteSpinor;
  typedef iImplPropagator<Simd>        SitePropagator;
@@ -711,6 +714,7 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
 
     typedef RealD  _Coeff_t ;
     static const int Dimension = Representation::Dimension;
+    static const bool isFundamental = Representation::isFundamental;
     static const bool LsVectorised=false;
     typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
       
@@ -839,6 +843,7 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
     public:
 
     static const int Dimension = Representation::Dimension;
+    static const bool isFundamental = Representation::isFundamental;
     static const bool LsVectorised=true;
     typedef RealD   Coeff_t ;
     typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
@@ -1033,29 +1038,29 @@ typedef WilsonImpl<vComplex,  TwoIndexAntiSymmetricRepresentation, CoeffReal > W
 typedef WilsonImpl<vComplexF, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplF;  // Float
 typedef WilsonImpl<vComplexD, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplD;  // Double
 
-typedef DomainWallVec5dImpl<vComplex ,Nc, CoeffReal> DomainWallVec5dImplR; // Real.. whichever prec
+typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation, CoeffReal> DomainWallVec5dImplR; // Real.. whichever prec
-typedef DomainWallVec5dImpl<vComplexF,Nc, CoeffReal> DomainWallVec5dImplF; // Float
+typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation, CoeffReal> DomainWallVec5dImplF; // Float
-typedef DomainWallVec5dImpl<vComplexD,Nc, CoeffReal> DomainWallVec5dImplD; // Double
+typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation, CoeffReal> DomainWallVec5dImplD; // Double
  
-typedef DomainWallVec5dImpl<vComplex ,Nc, CoeffRealHalfComms> DomainWallVec5dImplRL; // Real.. whichever prec
+typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation, CoeffRealHalfComms> DomainWallVec5dImplRL; // Real.. whichever prec
-typedef DomainWallVec5dImpl<vComplexF,Nc, CoeffRealHalfComms> DomainWallVec5dImplFH; // Float
+typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation, CoeffRealHalfComms> DomainWallVec5dImplFH; // Float
-typedef DomainWallVec5dImpl<vComplexD,Nc, CoeffRealHalfComms> DomainWallVec5dImplDF; // Double
+typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation, CoeffRealHalfComms> DomainWallVec5dImplDF; // Double
  
-typedef DomainWallVec5dImpl<vComplex ,Nc,CoeffComplex> ZDomainWallVec5dImplR; // Real.. whichever prec
+typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation,CoeffComplex> ZDomainWallVec5dImplR; // Real.. whichever prec
-typedef DomainWallVec5dImpl<vComplexF,Nc,CoeffComplex> ZDomainWallVec5dImplF; // Float
+typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation,CoeffComplex> ZDomainWallVec5dImplF; // Float
-typedef DomainWallVec5dImpl<vComplexD,Nc,CoeffComplex> ZDomainWallVec5dImplD; // Double
+typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation,CoeffComplex> ZDomainWallVec5dImplD; // Double
  
-typedef DomainWallVec5dImpl<vComplex ,Nc,CoeffComplexHalfComms> ZDomainWallVec5dImplRL; // Real.. whichever prec
+typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation,CoeffComplexHalfComms> ZDomainWallVec5dImplRL; // Real.. whichever prec
-typedef DomainWallVec5dImpl<vComplexF,Nc,CoeffComplexHalfComms> ZDomainWallVec5dImplFH; // Float
+typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation,CoeffComplexHalfComms> ZDomainWallVec5dImplFH; // Float
-typedef DomainWallVec5dImpl<vComplexD,Nc,CoeffComplexHalfComms> ZDomainWallVec5dImplDF; // Double
+typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation,CoeffComplexHalfComms> ZDomainWallVec5dImplDF; // Double
  
-typedef GparityWilsonImpl<vComplex , Nc,CoeffReal> GparityWilsonImplR;  // Real.. whichever prec
+typedef GparityWilsonImpl<vComplex , FundamentalRepresentation,CoeffReal> GparityWilsonImplR;  // Real.. whichever prec
-typedef GparityWilsonImpl<vComplexF, Nc,CoeffReal> GparityWilsonImplF;  // Float
+typedef GparityWilsonImpl<vComplexF, FundamentalRepresentation,CoeffReal> GparityWilsonImplF;  // Float
-typedef GparityWilsonImpl<vComplexD, Nc,CoeffReal> GparityWilsonImplD;  // Double
+typedef GparityWilsonImpl<vComplexD, FundamentalRepresentation,CoeffReal> GparityWilsonImplD;  // Double
  
-typedef GparityWilsonImpl<vComplex , Nc,CoeffRealHalfComms> GparityWilsonImplRL;  // Real.. whichever prec
+typedef GparityWilsonImpl<vComplex , FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplRL;  // Real.. whichever prec
-typedef GparityWilsonImpl<vComplexF, Nc,CoeffRealHalfComms> GparityWilsonImplFH;  // Float
+typedef GparityWilsonImpl<vComplexF, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplFH;  // Float
-typedef GparityWilsonImpl<vComplexD, Nc,CoeffRealHalfComms> GparityWilsonImplDF;  // Double
+typedef GparityWilsonImpl<vComplexD, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplDF;  // Double
 
 typedef StaggeredImpl<vComplex,  FundamentalRepresentation > StaggeredImplR;   // Real.. whichever prec
 typedef StaggeredImpl<vComplexF, FundamentalRepresentation > StaggeredImplF;  // Float

lib/qcd/action/fermion/WilsonFermion.cc

@@ -419,7 +419,7 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
     conformable(_grid, q_in._grid);
     conformable(_grid, q_out._grid);
     Lattice<iSinglet<Simd>> ph(_grid), coor(_grid);
-    Complex i(0.0,1.0);
+    ComplexD i(0.0,1.0);
     PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid);
     unsigned int tshift = (mu == Tp) ? 1 : 0;
     unsigned int LLt    = GridDefaultLatt()[Tp];
@@ -431,7 +431,7 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
         LatticeCoordinate(coor, mu);
         ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu])));
     }
-    ph = exp((Real)(2*M_PI)*i*ph);
+    ph = exp((RealD)(2*M_PI)*i*ph);
 
     q_out = zero;
     LatticeInteger coords(_grid);

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -793,7 +793,7 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
     Lattice<iSinglet<Simd>> ph(FermionGrid()), coor(FermionGrid());
     PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()),
                     tmp(FermionGrid());
-    Complex i(0.0, 1.0);
+    ComplexD i(0.0, 1.0);
     unsigned int tshift = (mu == Tp) ? 1 : 0;
     unsigned int LLs = q_in._grid->_rdimensions[0];
     unsigned int LLt    = GridDefaultLatt()[Tp];
@@ -806,7 +806,7 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
         LatticeCoordinate(coor, nu + 1);
         ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu])));
     }
-    ph = exp((Real)(2*M_PI)*i*ph);
+    ph = exp((RealD)(2*M_PI)*i*ph);
 
     q_out = zero;
     LatticeInteger coords(_FourDimGrid);

lib/qcd/action/fermion/WilsonKernels.h

@@ -55,7 +55,7 @@ template<class Impl> class WilsonKernels : public FermionOperator<Impl> , public
 public:
    
   template <bool EnableBool = true>
-  typename std::enable_if<Impl::Dimension == 3 && Nc == 3 &&EnableBool, void>::type
+  typename std::enable_if<Impl::isFundamental==true && Nc == 3 &&EnableBool, void>::type
   DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,
 		   int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) 
   {
@@ -99,7 +99,7 @@ public:
   }
      
   template <bool EnableBool = true>
-  typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type
+  typename std::enable_if<(Impl::isFundamental==false || (Impl::isFundamental==true && Nc != 3)) && EnableBool, void>::type
   DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,
 	   int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1 ) {
     // no kernel choice  
@@ -116,7 +116,7 @@ public:
   }
      
   template <bool EnableBool = true>
-  typename std::enable_if<Impl::Dimension == 3 && Nc == 3 && EnableBool,void>::type
+  typename std::enable_if<Impl::isFundamental==true && Nc == 3 && EnableBool,void>::type
   DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,
 	      int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) 
 {
@@ -161,7 +161,7 @@ public:
   }
 
   template <bool EnableBool = true>
-  typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool,void>::type
+  typename std::enable_if<(Impl::isFundamental==false || (Impl::isFundamental==true && Nc != 3)) && EnableBool,void>::type
   DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,SiteHalfSpinor * buf,
 		      int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) {
 

lib/qcd/action/scalar/Scalar.h

@@ -32,6 +32,7 @@ directory
 #include <Grid/qcd/action/scalar/ScalarImpl.h>
 #include <Grid/qcd/action/scalar/ScalarAction.h>
 #include <Grid/qcd/action/scalar/ScalarInteractionAction.h>
+#include <Grid/qcd/action/scalar/shGordonAction.h>
 
 namespace Grid {
 namespace QCD {
@@ -44,6 +45,9 @@ namespace QCD {
   template <int Colours, int Dimensions> using ScalarAdjActionF = ScalarInteractionAction<ScalarNxNAdjImplF<Colours>, Dimensions>;
   template <int Colours, int Dimensions> using ScalarAdjActionD = ScalarInteractionAction<ScalarNxNAdjImplD<Colours>, Dimensions>;
   
+  typedef shGordonAction<ScalarImplR>                 shGordonActionR;
+  typedef shGordonAction<ScalarImplF>                 shGordonActionF;
+  typedef shGordonAction<ScalarImplD>                 shGordonActionD;
 }
 }
 

lib/qcd/action/scalar/ScalarImpl.h

@@ -29,7 +29,9 @@ class ScalarImplTypes {
     static inline Field projectForce(Field& P){return P;}
 
     static inline void update_field(Field& P, Field& U, double ep) {
+      //std::cout << GridLogDebug << "P:\n" << P << std::endl;
       U += P*ep;
+      //std::cout << GridLogDebug << "U:\n" << U << std::endl;
     }
 
     static inline RealD FieldSquareNorm(Field& U) {
@@ -37,15 +39,17 @@ class ScalarImplTypes {
     }
 
     static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
-      gaussian(pRNG, U);
+     random(pRNG, U);
     }
 
     static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
-      gaussian(pRNG, U);
+      random(pRNG, U);
+      U *= 0.01;
     }
 
     static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
-      U = 1.0;
+      U = 0.0;
+      //std::cout << GridLogDebug << "Initial U:\n" << U << std::endl;
     }
 
     static void MomentumSpacePropagator(Field &out, RealD m)

lib/qcd/action/scalar/shGordonAction.h

@@ -0,0 +1,80 @@
+/*************************************************************************************
+
+  Grid physics library, www.github.com/paboyle/Grid
+
+  Source file: ./lib/qcd/action/gauge/shGordonAction.h
+
+  Copyright (C) 2018
+
+  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 */
+
+#ifndef SHGORDON_ACTION_H
+#define SHGORDON_ACTION_H
+
+namespace Grid {
+
+template <class Impl>
+class shGordonAction : public QCD::Action<typename Impl::Field> {
+ public:
+    INHERIT_FIELD_TYPES(Impl);
+
+ private:
+    RealD mass_square;
+    RealD g;
+
+ public:
+    shGordonAction(RealD ms, RealD g) : mass_square(ms), g(g) {}
+
+    virtual std::string LogParameters() {
+      std::stringstream sstream;
+      sstream << GridLogMessage << "[shGordonAction] g           : " << g           << std::endl;
+      sstream << GridLogMessage << "[shGordonAction] mass_square : " << mass_square << std::endl;
+      return sstream.str();
+    }
+    virtual std::string action_name() {return "shGordonAction";}
+
+    virtual void refresh(const Field &U, GridParallelRNG &pRNG) {}  // noop as no pseudoferms
+
+    virtual RealD S(const Field &phi) {
+      return QCD::Nd * ScalarObs<Impl>::sumphisquared(phi) + ScalarObs<Impl>::sumphider(phi) + 0.5*mass_square/(g*g)*sum(trace(exp(g*phi) + exp(-g*phi)))   ;
+    };
+
+    virtual void deriv(const Field &phi,
+                       Field &force) {
+        //std::cout << GridLogDebug << "Force total before :\n" << force << std::endl;
+        Field tmp(phi._grid);
+        tmp = 2.0*QCD::Nd*phi;
+        for (int mu = 0; mu < QCD::Nd; mu++) tmp -= Cshift(phi, mu, 1) + Cshift(phi, mu, -1);
+
+
+        //std::cout << GridLogDebug << "Phi norm : " << norm2(phi) << std::endl;
+        force += tmp + 0.5*mass_square/g*(exp(g*phi) - exp(-g*phi));
+        //std::cout << GridLogDebug << "Force tmp :\n" << tmp << std::endl;
+        //std::cout << GridLogDebug << "Force total after :\n" << force << std::endl;
+    }
+};
+
+
+
+}  // namespace Grid
+
+#endif // SHGORDON_ACTION_H

lib/qcd/hmc/integrators/Integrator.h

@@ -103,7 +103,7 @@ class Integrator {
         // Implement smearing only for the fundamental representation now
         repr_set.at(a)->deriv(Rep.U, forceR);
         GF force = Rep.RtoFundamentalProject(forceR);  // Ta for the fundamental rep
-        Real force_abs = std::sqrt(norm2(force)/(U._grid->gSites()));
+        Real force_abs = std::sqrt(norm2(force))/(U._grid->gSites());
         std::cout << GridLogIntegrator << "Hirep Force average: " << force_abs << std::endl;
         Mom -= force * ep ;
       }
@@ -116,6 +116,7 @@ class Integrator {
 
    for (int a = 0; a < as[level].actions.size(); ++a) {
       Field force(U._grid);
+      force = zero;
       conformable(U._grid, 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
@@ -123,7 +124,7 @@ class Integrator {
       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());
+      Real force_abs = std::sqrt(norm2(force))/U._grid->gSites();
       std::cout << GridLogIntegrator << "Force average: " << force_abs << std::endl;
       Mom -= force * ep; 
     }

lib/qcd/modules/ObservableModules.h

@@ -92,6 +92,20 @@ class PlaquetteMod: public ObservableModule<PlaquetteLogger<Impl>, NoParameters>
   PlaquetteMod(): ObsBase(NoParameters()){}
 };
 
+template < class Impl >
+class ExpScalarMod: public ObservableModule<ExpScalarLogger<Impl>, ExpScalarParameters>{
+  typedef ObservableModule<ExpScalarLogger<Impl>, ExpScalarParameters> ObsBase;
+  using ObsBase::ObsBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->ObservablePtr.reset(new ExpScalarLogger<Impl>(this->Par_));
+  }
+  public:
+  ExpScalarMod(ExpScalarParameters P): ObsBase(P){}
+  ExpScalarMod():ObsBase(){};
+};
+
 template < class Impl >
 class PolyakovMod: public ObservableModule<PolyakovLogger<Impl>, NoParameters>{
   typedef ObservableModule<PolyakovLogger<Impl>, NoParameters> ObsBase;

lib/qcd/observables/exp_scalar.h

@@ -0,0 +1,80 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/modules/exp_scalar.h
+
+Copyright (C) 2018
+
+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 */
+
+#ifndef HMC_EXP_SCALAR_H
+#define HMC_EXP_SCALAR_H
+
+namespace Grid {
+namespace QCD {
+
+struct ExpScalarParameters : Serializable {
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ExpScalarParameters,
+    double, a)
+
+    ExpScalarParameters(double _a = 0.0):a(_a){}
+
+    template < class ReaderClass >
+    ExpScalarParameters(Reader<ReaderClass>& Reader){
+        read(Reader, "ExpScalar", *this);  
+    }  
+};
+
+template <class Impl>
+class ExpScalarLogger : public HmcObservable<typename Impl::Field> {
+  ExpScalarParameters Pars;
+ public:
+
+  // necessary for HmcObservable compatibility
+  typedef typename Impl::Field Field;
+
+  ExpScalarLogger(double _a):Pars(_a){}
+
+  ExpScalarLogger(ExpScalarParameters P):Pars(P){}
+
+  void TrajectoryComplete(int traj, typename Impl::Field &U,
+                          GridSerialRNG &sRNG,
+                          GridParallelRNG &pRNG) {
+
+    double e = sum(trace(exp(Pars.a*U)));
+
+    int def_prec = std::cout.precision();
+
+    std::cout << GridLogMessage
+        << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
+        << "ExpScalar: [ " << traj << " ] "<< e << std::endl;
+
+    std::cout.precision(def_prec);
+
+  }
+};
+
+}  // namespace QCD
+}  // namespace Grid
+
+#endif  // HMC_PLAQUETTE_H

lib/qcd/observables/hmc_observable.h

@@ -46,6 +46,6 @@ class HmcObservable {
 #include "plaquette.h"
 #include "topological_charge.h"
 #include "polyakov_loop.h"
-
+#include "exp_scalar.h"
 
 #endif  //  HMC_OBSERVABLE_H

lib/qcd/representations/adjoint.h

@@ -23,6 +23,7 @@ class AdjointRep {
   typedef typename SU_Adjoint<ncolour>::LatticeAdjMatrix LatticeMatrix;
   typedef typename SU_Adjoint<ncolour>::LatticeAdjField LatticeField;
   static const int Dimension = ncolour * ncolour - 1;
+  static const bool isFundamental = false;
 
   LatticeField U;
 

lib/qcd/representations/fundamental.h

@@ -19,6 +19,7 @@ template <int ncolour>
 class FundamentalRep {
  public:
   static const int Dimension = ncolour;
+  static const bool isFundamental = true;
 
   // typdef to be used by the Representations class in HMC to get the
   // types for the higher representation fields

lib/qcd/representations/two_index.h

@@ -29,6 +29,7 @@ class TwoIndexRep {
   typedef typename SU_TwoIndex<ncolour, S>::LatticeTwoIndexMatrix LatticeMatrix;
   typedef typename SU_TwoIndex<ncolour, S>::LatticeTwoIndexField LatticeField;
   static const int Dimension = ncolour * (ncolour + S) / 2;
+  static const bool isFundamental = false;
 
   LatticeField U;
 

lib/version.h

@@ -0,0 +1 @@
+

tests/hmc/Test_hmc_shGordonAction.cc

@@ -0,0 +1,118 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_shGordonAction.cc
+
+Copyright (C) 2018
+
+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 {
+class ScalarActionParameters : Serializable {
+ public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarActionParameters,
+    double, mass_squared,
+    double, g);
+
+    template <class ReaderClass >
+  ScalarActionParameters(Reader<ReaderClass>& Reader){
+    read(Reader, "ScalarAction", *this);
+  }
+
+};
+}
+
+
+
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  typedef Grid::JSONReader       Serialiser;
+
+  Grid_init(&argc, &argv);
+  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 ScalarGenericHMCRunner HMCWrapper;  // Uses the default minimum norm, real scalar fields
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+  TheHMC.ReadCommandLine(argc, argv);
+
+  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);
+  // Grid from the command line
+  constexpr int Ndimensions = 2;
+  GridModule ScalarGrid;
+  if (GridDefaultLatt().size() != Ndimensions){
+    std::cout << "Incorrect dimension of the grid\n. Expected dim="<< Ndimensions << std::endl;
+    exit(1);
+  }
+  if (GridDefaultMpi().size() != Ndimensions){
+    std::cout << "Incorrect dimension of the mpi grid\n. Expected dim="<< Ndimensions << std::endl;
+    exit(1);
+  }
+  ScalarGrid.set_full(new GridCartesian(GridDefaultLatt(),GridDefaultSimd(Ndimensions, vComplex::Nsimd()),GridDefaultMpi()));
+  ScalarGrid.set_rb(new GridRedBlackCartesian(ScalarGrid.get_full()));
+  TheHMC.Resources.AddGrid("scalar", ScalarGrid);
+  std::cout << "Lattice size : " << GridDefaultLatt() << std::endl;
+
+  CheckpointerParameters CPparams(Reader);
+  TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar(Reader);
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Some online observable measurements
+  typedef ExpScalarMod<HMCWrapper::ImplPolicy> ExpObs;
+  ExpScalarParameters ExpParams(Reader);
+  TheHMC.Resources.AddObservable<ExpObs>(ExpParams);
+  ///////////////////////////////////////////
+
+  // Real Scalar sh-Gordon action
+  ScalarActionParameters SPar(Reader);
+  shGordonActionR Saction(SPar.mass_squared, SPar.g);
+
+  // Collect actions
+  ActionLevel<shGordonActionR::Field, ScalarFields> Level1(1);
+  Level1.push_back(&Saction);
+
+  TheHMC.TheAction.push_back(Level1);
+  /////////////////////////////////////////////////////////////
+
+  TheHMC.Parameters.initialize(Reader);
+  TheHMC.Run(); 
+
+  Grid_finalize();
+
+} // main
+
+