Fix for the CI

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/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
+
+