First implementation of the sh-Gordon action

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/shGordonAction.h

@@ -0,0 +1,75 @@
+/*************************************************************************************
+
+  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 0.5*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) {
+        Field tmp(phi._grid);
+        tmp = zero;
+        for (int mu = 0; mu < QCD::Nd; mu++) tmp += Cshift(phi, mu, 1) - Cshift(phi, mu, -1);
+
+        force+= 0.5*tmp - 0.5*mass_square/g*(exp(g*phi) - exp(-g*phi));
+    }
+};
+
+
+
+}  // namespace Grid
+
+#endif // SHGORDON_ACTION_H

tests/hmc/Test_hmc_shGordonAction.cc

@@ -0,0 +1,100 @@
+/*************************************************************************************
+
+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>
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  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
+  //typedef Representations<EmptyRep<typename ScalarMatrixImplTypes<vComplex, 3>::Field> > ScalarMatrixFields;
+  //typedef HMCWrapperTemplate<ScalarMatrixImplTypes<vComplex, 3>, MinimumNorm2, ScalarMatrixFields> HMCWrapper;
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+
+  // Grid from the command line
+  GridModule ScalarGrid;
+  ScalarGrid.set_full( SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(), GridDefaultSimd(Nd, vReal::Nsimd()),
+        GridDefaultMpi()));
+  ScalarGrid.set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(ScalarGrid.get_full()));
+  TheHMC.Resources.AddGrid("scalar", ScalarGrid);
+  // Possibile to create the module by hand 
+  // hardcoding parameters or using a Reader
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams;  
+  CPparams.config_prefix = "ckpoint_scalar_lat";
+  CPparams.rng_prefix = "ckpoint_scalar_rng";
+  CPparams.saveInterval = 50;
+  CPparams.format = "IEEE64BIG";
+  
+  TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar;
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes 
+  // that have a complex construction
+  // standard
+
+  // Real Scalar action
+  shGordonActionR Saction(1.0,1.0);
+
+  // Collect actions
+  ActionLevel<shGordonActionR::Field, ScalarFields> Level1(1);
+  Level1.push_back(&Saction);
+
+  TheHMC.TheAction.push_back(Level1);
+  /////////////////////////////////////////////////////////////
+
+
+  // HMC parameters are serialisable 
+  TheHMC.Parameters.MD.MDsteps = 60;
+  TheHMC.Parameters.MD.trajL   = 1.0;
+
+  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  TheHMC.Run(); 
+
+  Grid_finalize();
+
+} // main
+
+