From e6acffdfc21c5e02ac3836abc43ad7324aa534b0 Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Fri, 21 Oct 2016 16:06:34 +0100
Subject: [PATCH] Fixing the plaquette computation

---
 lib/qcd/action/gauge/WilsonGaugeAction.h | 11 +++-------
 lib/qcd/utils/SUn.h                      |  3 +++
 lib/qcd/utils/WilsonLoops.h              | 12 ++++++++++
 tests/hmc/Test_hmc_WilsonGauge_Binary.cc | 28 ++++++++++--------------
 4 files changed, 29 insertions(+), 25 deletions(-)

diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h
index 059fcc58..8ac1df74 100644
--- a/lib/qcd/action/gauge/WilsonGaugeAction.h
+++ b/lib/qcd/action/gauge/WilsonGaugeAction.h
@@ -99,7 +99,7 @@ class VariableWilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
 
  public:
   VariableWilsonGaugeAction(std::vector<RealD> bulk, std::vector<RealD> xdim,
-                            GridBase *_grid, bool openBC = false)
+                            GridBase *_grid)
       : b_bulk(bulk),
         b_xdim(xdim),
         grid(_grid),
@@ -130,11 +130,6 @@ class VariableWilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
       beta_xdim = where(coor == tau, temp, beta_xdim);
     }
 
-    if (!openBC) {
-      temp = b_xdim[Nex - 1];
-      beta_xdim = where(coor == Nex - 1, temp, beta_xdim);
-    }
-
     beta_xdim_shifted = Cshift(beta_xdim, Ndim - 1, -1);
 
     beta_bulk = zero;
@@ -143,8 +138,7 @@ class VariableWilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
       beta_bulk = where(coor == tau, temp, beta_bulk);
     }
 
-    std::cout << beta_xdim << std::endl;
-    std::cout << beta_xdim_shifted << std::endl;
+
   };
 
   virtual void refresh(const GaugeField &U,
@@ -180,6 +174,7 @@ class VariableWilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
       }
     }
 
+    
     double faces = (1.0 * (Nd - 1) * (Nd - 2)) / 2.0;
     SumdirPlaq *= OneOnNc / (RealD(bulk_volume) * faces);
 
diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h
index 20cd1889..9ad3ca18 100644
--- a/lib/qcd/utils/SUn.h
+++ b/lib/qcd/utils/SUn.h
@@ -170,6 +170,7 @@ class SU {
     ta()()(i2, i1) = 1.0;
     ta = ta * 0.5;
   }
+  
   template <class cplx>
   static void generatorSigmaX(int su2Index, iSUnMatrix<cplx> &ta) {
     ta = zero;
@@ -194,6 +195,8 @@ class SU {
     ta = ta * nrm;
   }
 
+
+
   ////////////////////////////////////////////////////////////////////////
   // Map a su2 subgroup number to the pair of rows that are non zero
   ////////////////////////////////////////////////////////////////////////
diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h
index e47040ee..19de7ab0 100644
--- a/lib/qcd/utils/WilsonLoops.h
+++ b/lib/qcd/utils/WilsonLoops.h
@@ -54,9 +54,21 @@ public:
     // resolution throughout the usage in this file, and rather defeats the
     // purpose of deriving
     // from Gimpl.
+    /*
     plaq = Gimpl::CovShiftBackward(
         U[mu], mu, Gimpl::CovShiftBackward(
                        U[nu], nu, Gimpl::CovShiftForward(U[mu], mu, U[nu])));
+                       */
+    // _
+    //|< _|
+    plaq = Gimpl::CovShiftForward(U[mu],mu,
+           Gimpl::CovShiftForward(U[nu],nu,
+           Gimpl::CovShiftBackward(U[mu],mu,
+           Gimpl::CovShiftIdentityBackward(U[nu], nu))));
+
+
+
+
   }
   //////////////////////////////////////////////////
   // trace of directed plaquette oriented in mu,nu plane
diff --git a/tests/hmc/Test_hmc_WilsonGauge_Binary.cc b/tests/hmc/Test_hmc_WilsonGauge_Binary.cc
index 9a9c3e8e..a3bdbe4f 100644
--- a/tests/hmc/Test_hmc_WilsonGauge_Binary.cc
+++ b/tests/hmc/Test_hmc_WilsonGauge_Binary.cc
@@ -62,32 +62,26 @@ class HmcRunner : public BinaryHmcRunner {
   void BuildTheAction(int argc, char **argv)
 
   {
+    int Ndim=5;
     typedef WilsonImplR ImplPolicy;
     typedef WilsonFermionR FermionAction;
     typedef typename FermionAction::FermionField FermionField;
 
-    UGrid = SpaceTimeGrid::makeFourDimGrid(
-        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
-        GridDefaultMpi());
-    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+    std::vector<int> simd = GridDefaultSimd(Ndim-1,vComplex::Nsimd());
+    simd.push_back(1);
 
-    FGrid = UGrid;
-    FrbGrid = UrbGrid;
 
-    // temporarily need a gauge field
-    LatticeGaugeField U(UGrid);
+    //UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), simd, GridDefaultMpi());
+    UGrid   =  new GridCartesian(GridDefaultLatt(),simd,GridDefaultMpi());
 
     // Gauge action
-    int Ls = UGrid->_gdimensions[Nd - 1];
-    std::vector<RealD> betat(Ls);
-    std::vector<RealD> betas(Ls);
-    //betat={5,6,6,6,6,6,6,5};
-    betat={1,1,1,0,1,1,1,1};
-    //betas={5.2,5.5,5.8,6,6,5.8,5.5,5.2};
-    betas={0,0,0,0,0,0,0,0};
-    bool openBC = false;
+    int Ls = UGrid->_fdimensions[Nd - 1];
+    std::vector<RealD> betat(Ls,6.0);
+    std::vector<RealD> betas(Ls,5.6);
+    betat[Ls-1]= 0.0;
+    betas={5.2,5.5,5.8,6,6,5.8,5.5,5.2};
     std:cout << GridLogMessage << "Betas: " << betas << std::endl;
-    VariableWilsonGaugeActionR Waction(betas, betat, UGrid, openBC);
+    VariableWilsonGaugeActionR Waction(betas, betat, UGrid);
     //WilsonGaugeActionR Waction(5.6);
 
     // Collect actions