From f5cc7e253bf21eda17f8c74cd850ba30c2b34bb3 Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Thu, 18 Jan 2018 10:24:28 +0000
Subject: [PATCH] Debugging the shGordon action

---
 lib/qcd/QCD.h                          |  2 +-
 lib/qcd/action/scalar/shGordonAction.h |  8 ++++----
 tests/hmc/Test_hmc_shGordonAction.cc   | 21 +++++++++++++++++++--
 3 files changed, 24 insertions(+), 7 deletions(-)

diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h
index e1e38e82..16f82af1 100644
--- a/lib/qcd/QCD.h
+++ b/lib/qcd/QCD.h
@@ -51,7 +51,7 @@ namespace QCD {
 
     static const int Nc=3;
     static const int Ns=4;
-    static const int Nd=4;
+    static const int Nd=2;
     static const int Nhs=2; // half spinor
     static const int Nds=8; // double stored gauge field
     static const int Ngp=2; // gparity index range
diff --git a/lib/qcd/action/scalar/shGordonAction.h b/lib/qcd/action/scalar/shGordonAction.h
index 3b402c72..12f208ac 100644
--- a/lib/qcd/action/scalar/shGordonAction.h
+++ b/lib/qcd/action/scalar/shGordonAction.h
@@ -55,16 +55,16 @@ class shGordonAction : public QCD::Action<typename Impl::Field> {
     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)))   ;
+      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) {
         Field tmp(phi._grid);
-        tmp = zero;
-        for (int mu = 0; mu < QCD::Nd; mu++) tmp += Cshift(phi, mu, 1) - Cshift(phi, mu, -1);
+        tmp = 2.0*QCD::Nd*phi;
+        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));
+        force+= tmp - 0.5*mass_square/g*(exp(g*phi) - exp(-g*phi));
     }
 };
 
diff --git a/tests/hmc/Test_hmc_shGordonAction.cc b/tests/hmc/Test_hmc_shGordonAction.cc
index 96be6a96..45cb6b62 100644
--- a/tests/hmc/Test_hmc_shGordonAction.cc
+++ b/tests/hmc/Test_hmc_shGordonAction.cc
@@ -44,12 +44,29 @@ int main(int argc, char **argv) {
   HMCWrapper TheHMC;
 
   // 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;
+
+  /*
   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
 
@@ -76,7 +93,7 @@ int main(int argc, char **argv) {
   // standard
 
   // Real Scalar action
-  shGordonActionR Saction(1.0,1.0);
+  shGordonActionR Saction(0.1,0.1);
 
   // Collect actions
   ActionLevel<shGordonActionR::Field, ScalarFields> Level1(1);
@@ -87,7 +104,7 @@ int main(int argc, char **argv) {
 
 
   // HMC parameters are serialisable 
-  TheHMC.Parameters.MD.MDsteps = 60;
+  TheHMC.Parameters.MD.MDsteps = 100;
   TheHMC.Parameters.MD.trajL   = 1.0;
 
   TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file