From c2250fa124e90a2fc255e3ff043ab3935364a8a7 Mon Sep 17 00:00:00 2001
From: Nils Asmussen <n.asmussen@soton.ac.uk>
Date: Fri, 29 Mar 2019 16:36:56 +0000
Subject: [PATCH 1/8] MFermion::GaugeProp fix for 4d fields

---
 Hadrons/Modules/MFermion/GaugeProp.hpp | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/Hadrons/Modules/MFermion/GaugeProp.hpp b/Hadrons/Modules/MFermion/GaugeProp.hpp
index f11113b4..49b4bd5a 100644
--- a/Hadrons/Modules/MFermion/GaugeProp.hpp
+++ b/Hadrons/Modules/MFermion/GaugeProp.hpp
@@ -111,13 +111,18 @@ void TGaugeProp<FImpl>::setup(void)
 {
     Ls_ = env().getObjectLs(par().solver);
     envCreateLat(PropagatorField, getName());
-    envTmpLat(FermionField, "source", Ls_);
-    envTmpLat(FermionField, "sol", Ls_);
     envTmpLat(FermionField, "tmp");
     if (Ls_ > 1)
     {
+        envTmpLat(FermionField, "source", Ls_);
+        envTmpLat(FermionField, "sol", Ls_);
         envCreateLat(PropagatorField, getName() + "_5d", Ls_);
     }
+    else
+    {
+       envTmpLat(FermionField, "source");
+       envTmpLat(FermionField, "sol");
+    }
 }
 
 // execution ///////////////////////////////////////////////////////////////////

From edeb5908180088cf3e17280bfffd4bbbbabca2a6 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Fri, 3 May 2019 17:09:47 +0100
Subject: [PATCH 2/8] DiskVector: fix of memory bug triggering segfault when
 the cache is accessed following a certain pattern

---
 Hadrons/DiskVector.hpp | 26 ++++++++++++++++++++------
 1 file changed, 20 insertions(+), 6 deletions(-)

diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp
index 7410f1c3..3a8d0696 100644
--- a/Hadrons/DiskVector.hpp
+++ b/Hadrons/DiskVector.hpp
@@ -395,12 +395,26 @@ void DiskVectorBase<T>::cacheInsert(const unsigned int i, const T &obj) const
     auto &freeInd  = *freePtr_;
     auto &loads    = *loadsPtr_;
 
-    evict();
-    index[i] = freeInd.top();
-    freeInd.pop();
-    cache[index.at(i)] = obj;
-    loads.push_back(i);
-    modified[index.at(i)] = false;
+    // cache miss, evict and store
+    if (index.find(i) == index.end())
+    {
+        evict();
+        index[i] = freeInd.top();
+        freeInd.pop();
+        cache[index.at(i)] = obj;
+        loads.push_back(i);
+        modified[index.at(i)] = false;
+    }
+    // cache hit, modify current value
+    else
+    {
+        auto pos = std::find(loads.begin(), loads.end(), i);
+        
+        cache[index.at(i)]    = obj;
+        modified[index.at(i)] = true;
+        loads.erase(pos);
+        loads.push_back(i);
+    }
 
 #ifdef DV_DEBUG
     std::string msg;

From 2acd8ece655af8882195c6581c16712d21cc65fc Mon Sep 17 00:00:00 2001
From: fionnoh <fionnoh@gmail.com>
Date: Wed, 8 May 2019 10:57:36 +0100
Subject: [PATCH 3/8] Hadron WeakEye and A2ALoop bug fixes, and WWVVContraction
 bug fix

---
 Grid/qcd/utils/A2Autils.h                   | 9 +++++----
 Hadrons/Modules/MContraction/A2ALoop.hpp    | 2 +-
 Hadrons/Modules/MContraction/WeakEye3pt.hpp | 6 +++---
 3 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h
index 89b4d4bd..97188ffe 100644
--- a/Grid/qcd/utils/A2Autils.h
+++ b/Grid/qcd/utils/A2Autils.h
@@ -986,17 +986,18 @@ void A2Autils<FImpl>::ContractWWVV(std::vector<PropagatorField> &WWVV,
       for(int t=0;t<N_t;t++){
       for(int s=0;s<N_s;s++){
 	auto tmp1 = vs[s]._odata[ss];
-	vobj tmp2 = zero;
+  vobj tmp2 = zero;
+  vobj tmp3 = zero;
 
 	for(int d=d_o;d<MIN(d_o+d_unroll,N_d);d++){
 	  Scalar_v coeff = WW_sd(t,s,d);
-	  mac(&tmp2 ,& coeff, & vd[d]._odata[ss]);
-	}
+	  tmp3 = conjugate(vd[d]._odata[ss]);
+	  mac(&tmp2, &coeff, &tmp3);
+  }
 
 	//////////////////////////
 	// Fast outer product of tmp1 with a sum of terms suppressed by d_unroll
 	//////////////////////////
-	tmp2 = conjugate(tmp2);
 	for(int s1=0;s1<Ns;s1++){
 	for(int s2=0;s2<Ns;s2++){
 	  WWVV[t]._odata[ss]()(s1,s2)(0,0) += tmp1()(s1)(0)*tmp2()(s2)(0);
diff --git a/Hadrons/Modules/MContraction/A2ALoop.hpp b/Hadrons/Modules/MContraction/A2ALoop.hpp
index 8ae129e3..2ef99354 100644
--- a/Hadrons/Modules/MContraction/A2ALoop.hpp
+++ b/Hadrons/Modules/MContraction/A2ALoop.hpp
@@ -89,7 +89,7 @@ std::vector<std::string> TA2ALoop<FImpl>::getInput(void)
 template <typename FImpl>
 std::vector<std::string> TA2ALoop<FImpl>::getOutput(void)
 {
-    std::vector<std::string> out = {};
+    std::vector<std::string> out = {getName()};
     
     return out;
 }
diff --git a/Hadrons/Modules/MContraction/WeakEye3pt.hpp b/Hadrons/Modules/MContraction/WeakEye3pt.hpp
index 22af2501..ea7ff529 100644
--- a/Hadrons/Modules/MContraction/WeakEye3pt.hpp
+++ b/Hadrons/Modules/MContraction/WeakEye3pt.hpp
@@ -52,7 +52,7 @@ BEGIN_HADRONS_NAMESPACE
  *                                        |
  *                one trace               |                two traces
  * 
- * one trace : tr(qbr*gOut*qs*adj(gIn)*g5*adj(qbl)*g5*G*loop*G*qbr*gOut)
+ * one trace : tr(qbr*gOut*qs*adj(gIn)*g5*adj(qbl)*g5*G*loop*G)
  * two traces: tr(qbr*gOut*qs*adj(gIn)*g5*adj(qbl)*g5*G)*tr(loop*G)
  * 
  */
@@ -118,7 +118,7 @@ template <typename FImpl>
 std::vector<std::string> TWeakEye3pt<FImpl>::getInput(void)
 {
     std::vector<std::string> in = {par().qBarLeft, par().qBarRight, 
-                                   par().qSpectator};
+                                   par().qSpectator, par().loop};
     
     return in;
 }
@@ -170,7 +170,7 @@ void TWeakEye3pt<FImpl>::execute(void)
 
         r.info.op = G.g;
         // one trace
-        corr = trace(qbr*gOut*qst*adj(gIn)*g5*adj(qbl)*g5*G*loop*G*qbr*gOut);
+        corr = trace(qbr*gOut*qst*adj(gIn)*g5*adj(qbl)*g5*G*loop*G);
         sliceSum(corr, buf, Tp);
         r.corr.clear();
         for (unsigned int t = 0; t < buf.size(); ++t)

From 94accec311acfa99f6ce7304b7fcbc6476dc3e69 Mon Sep 17 00:00:00 2001
From: fionnoh <fionnoh@pscomp.ie>
Date: Wed, 15 May 2019 13:35:47 +0100
Subject: [PATCH 4/8] Added gauge transform option to eigpack IO

---
 Hadrons/Modules/MIO/LoadEigenPack.cc  |  4 +-
 Hadrons/Modules/MIO/LoadEigenPack.hpp | 64 +++++++++++++++++++++------
 2 files changed, 52 insertions(+), 16 deletions(-)

diff --git a/Hadrons/Modules/MIO/LoadEigenPack.cc b/Hadrons/Modules/MIO/LoadEigenPack.cc
index 15153871..0750d265 100644
--- a/Hadrons/Modules/MIO/LoadEigenPack.cc
+++ b/Hadrons/Modules/MIO/LoadEigenPack.cc
@@ -31,7 +31,7 @@ using namespace Grid;
 using namespace Hadrons;
 using namespace MIO;
 
-template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL>>;
+template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL>, GIMPL>;
 #ifdef GRID_DEFAULT_PRECISION_DOUBLE
-template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>>;
+template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>, GIMPL>;
 #endif
diff --git a/Hadrons/Modules/MIO/LoadEigenPack.hpp b/Hadrons/Modules/MIO/LoadEigenPack.hpp
index 9751073a..1fffb645 100644
--- a/Hadrons/Modules/MIO/LoadEigenPack.hpp
+++ b/Hadrons/Modules/MIO/LoadEigenPack.hpp
@@ -47,16 +47,21 @@ public:
                                     std::string, filestem,
                                     bool, multiFile,
                                     unsigned int, size,
-                                    unsigned int, Ls);
+                                    unsigned int, Ls,
+                                    std::string, gaugeXform);
 };
 
-template <typename Pack>
+template <typename Pack, typename GImpl>
 class TLoadEigenPack: public Module<LoadEigenPackPar>
 {
 public:
     typedef typename Pack::Field   Field;
     typedef typename Pack::FieldIo FieldIo;
     typedef BaseEigenPack<Field>   BasePack;
+
+public:
+    GAUGE_TYPE_ALIASES(GImpl, );
+    typedef typename GImpl::GaugeLinkField GaugeMat;
 public:
     // constructor
     TLoadEigenPack(const std::string name);
@@ -71,31 +76,36 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_TMP(LoadFermionEigenPack, TLoadEigenPack<FermionEigenPack<FIMPL>>, MIO);
+MODULE_REGISTER_TMP(LoadFermionEigenPack, ARG(TLoadEigenPack<FermionEigenPack<FIMPL>, GIMPL>), MIO);
 #ifdef GRID_DEFAULT_PRECISION_DOUBLE
-MODULE_REGISTER_TMP(LoadFermionEigenPackIo32, ARG(TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>>), MIO);
+MODULE_REGISTER_TMP(LoadFermionEigenPackIo32, ARG(TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>, GIMPL>), MIO);
 #endif
 
 /******************************************************************************
  *                    TLoadEigenPack implementation                           *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-template <typename Pack>
-TLoadEigenPack<Pack>::TLoadEigenPack(const std::string name)
+template <typename Pack, typename GImpl>
+TLoadEigenPack<Pack, GImpl>::TLoadEigenPack(const std::string name)
 : Module<LoadEigenPackPar>(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
-template <typename Pack>
-std::vector<std::string> TLoadEigenPack<Pack>::getInput(void)
+template <typename Pack, typename GImpl>
+std::vector<std::string> TLoadEigenPack<Pack, GImpl>::getInput(void)
 {
     std::vector<std::string> in;
+
+    if (!par().gaugeXform.empty())
+    {
+        in = {par().gaugeXform};
+    }
     
     return in;
 }
 
-template <typename Pack>
-std::vector<std::string> TLoadEigenPack<Pack>::getOutput(void)
+template <typename Pack, typename GImpl>
+std::vector<std::string> TLoadEigenPack<Pack, GImpl>::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
     
@@ -103,8 +113,8 @@ std::vector<std::string> TLoadEigenPack<Pack>::getOutput(void)
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
-template <typename Pack>
-void TLoadEigenPack<Pack>::setup(void)
+template <typename Pack, typename GImpl>
+void TLoadEigenPack<Pack, GImpl>::setup(void)
 {
     GridBase *gridIo = nullptr;
 
@@ -114,16 +124,42 @@ void TLoadEigenPack<Pack>::setup(void)
     }
     envCreateDerived(BasePack, Pack, getName(), par().Ls, par().size, 
                      envGetRbGrid(Field, par().Ls), gridIo);
+
+    if (!par().gaugeXform.empty())
+    {
+        envTmp(GaugeMat,    "tmp5dXform", par().Ls, envGetGrid5(Field, par().Ls));
+        envTmp(GaugeMat, "tmp5dXformOdd", par().Ls, envGetRbGrid5(Field, par().Ls));
+    }
 }
 
 // execution ///////////////////////////////////////////////////////////////////
-template <typename Pack>
-void TLoadEigenPack<Pack>::execute(void)
+template <typename Pack, typename GImpl>
+void TLoadEigenPack<Pack, GImpl>::execute(void)
 {
     auto &epack = envGetDerived(BasePack, Pack, getName());
 
     epack.read(par().filestem, par().multiFile, vm().getTrajectory());
     epack.eval.resize(par().size);
+
+    if (!par().gaugeXform.empty())
+    {
+        auto &xform = envGet(GaugeMat, par().gaugeXform);
+        envGetTmp(GaugeMat,    tmp5dXform);
+        envGetTmp(GaugeMat, tmp5dXformOdd);
+
+        for (unsigned int j = 0; j < par().Ls; j++)
+        {
+            InsertSlice(xform, tmp5dXform, j, 0);
+        }
+
+        pickCheckerboard(Odd, tmp5dXformOdd, tmp5dXform);
+        for (unsigned int i = 0; i < par().size; i++)
+        {
+            LOG(Message) << "Applying gauge transformation to eigenvector i = " << i << std::endl;
+            epack.evec[i].checkerboard = Odd;
+            epack.evec[i] = tmp5dXformOdd * epack.evec[i];
+        }
+    }
 }
 
 END_MODULE_NAMESPACE

From ce102ac550f0d2ca28f408b06cb86c606d5f3a19 Mon Sep 17 00:00:00 2001
From: fionnoh <fionnoh@pscomp.ie>
Date: Wed, 15 May 2019 14:31:25 +0100
Subject: [PATCH 5/8] More logging, timing, and 4d/5d logic for eigpack gauge
 transforms

---
 Hadrons/Modules/MIO/LoadEigenPack.hpp | 42 ++++++++++++++++++++-------
 1 file changed, 32 insertions(+), 10 deletions(-)

diff --git a/Hadrons/Modules/MIO/LoadEigenPack.hpp b/Hadrons/Modules/MIO/LoadEigenPack.hpp
index 1fffb645..0d6036b8 100644
--- a/Hadrons/Modules/MIO/LoadEigenPack.hpp
+++ b/Hadrons/Modules/MIO/LoadEigenPack.hpp
@@ -127,8 +127,19 @@ void TLoadEigenPack<Pack, GImpl>::setup(void)
 
     if (!par().gaugeXform.empty())
     {
-        envTmp(GaugeMat,    "tmp5dXform", par().Ls, envGetGrid5(Field, par().Ls));
-        envTmp(GaugeMat, "tmp5dXformOdd", par().Ls, envGetRbGrid5(Field, par().Ls));
+        if (par().Ls > 1)
+        {
+            LOG(Message) << "Setup 5d GaugeMat for Ls = " << par().Ls << std::endl;
+            envTmp(GaugeMat,    "tmpXform", par().Ls, envGetGrid5(Field, par().Ls));
+            envTmp(GaugeMat, "tmpXformOdd", par().Ls, envGetRbGrid5(Field, par().Ls));
+        }
+        else
+        {
+            LOG(Message) << "Setup 4d GaugeMat for Ls = " << par().Ls << std::endl;
+            envTmp(GaugeMat,    "tmpXform", par().Ls, envGetGrid(Field));
+            envTmp(GaugeMat, "tmpXformOdd", par().Ls, envGetRbGrid(Field));
+        }
+        
     }
 }
 
@@ -143,22 +154,33 @@ void TLoadEigenPack<Pack, GImpl>::execute(void)
 
     if (!par().gaugeXform.empty())
     {
-        auto &xform = envGet(GaugeMat, par().gaugeXform);
-        envGetTmp(GaugeMat,    tmp5dXform);
-        envGetTmp(GaugeMat, tmp5dXformOdd);
 
-        for (unsigned int j = 0; j < par().Ls; j++)
+        LOG(Message) << "Applying gauge transformation to eigenvectors " << getName()
+                     << " using " << par().gaugeXform << std::endl;
+        auto &xform = envGet(GaugeMat, par().gaugeXform);
+        envGetTmp(GaugeMat,    tmpXform);
+        envGetTmp(GaugeMat, tmpXformOdd);
+
+        if (par().Ls > 1) 
         {
-            InsertSlice(xform, tmp5dXform, j, 0);
+            LOG(Message) << "Creating 5d GaugeMat from " << par().gaugeXform << std::endl;
+            startTimer("5-d gauge transform creation");
+            for (unsigned int j = 0; j < par().Ls; j++)
+            {
+                InsertSlice(xform, tmpXform, j, 0);
+            }
+            stopTimer("5-d gauge transform creation");
         }
 
-        pickCheckerboard(Odd, tmp5dXformOdd, tmp5dXform);
+        pickCheckerboard(Odd, tmpXformOdd, tmpXform);
+        startTimer("Transform application");
         for (unsigned int i = 0; i < par().size; i++)
         {
-            LOG(Message) << "Applying gauge transformation to eigenvector i = " << i << std::endl;
+            LOG(Message) << "Applying gauge transformation to eigenvector i = " << i << "/" << par().size << std::endl;
             epack.evec[i].checkerboard = Odd;
-            epack.evec[i] = tmp5dXformOdd * epack.evec[i];
+            epack.evec[i] = tmpXformOdd * epack.evec[i];
         }
+        stopTimer("Transform application");
     }
 }
 

From 48b1c806edbdab3294f204ec582ee0e0d11db6b2 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Fri, 17 May 2019 17:36:32 +0100
Subject: [PATCH 6/8] Coulomb gauge added as an option

---
 Grid/qcd/utils/GaugeFix.h   | 72 +++++++++++++++++++++++++------------
 tests/core/Test_fft_gfix.cc | 20 +++++++++++
 2 files changed, 70 insertions(+), 22 deletions(-)

diff --git a/Grid/qcd/utils/GaugeFix.h b/Grid/qcd/utils/GaugeFix.h
index 2b0cd7f2..67f84f14 100644
--- a/Grid/qcd/utils/GaugeFix.h
+++ b/Grid/qcd/utils/GaugeFix.h
@@ -31,6 +31,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 namespace Grid {
 namespace QCD {
 
+
 template <class Gimpl> 
 class FourierAcceleratedGaugeFixer  : public Gimpl {
  public:
@@ -45,18 +46,21 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
       A[mu] = Ta(U[mu]) * cmi;
     }
   }
-  static void DmuAmu(const std::vector<GaugeMat> &A,GaugeMat &dmuAmu) {
+  static void DmuAmu(const std::vector<GaugeMat> &A,GaugeMat &dmuAmu,int orthog) {
     dmuAmu=zero;
     for(int mu=0;mu<Nd;mu++){
-      dmuAmu = dmuAmu + A[mu] - Cshift(A[mu],mu,-1);
+      if ( mu != orthog ) {
+	dmuAmu = dmuAmu + A[mu] - Cshift(A[mu],mu,-1);
+      }
     }
   }  
-  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false) {
+
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
     GridBase *grid = Umu._grid;
     GaugeMat xform(grid);
-    SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier);
+    SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog);
   }
-  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false) {
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
 
     GridBase *grid = Umu._grid;
 
@@ -71,15 +75,29 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
 
     GaugeMat dmuAmu(grid);
 
+    {
+      Real plaq      =WilsonLoops<Gimpl>::avgPlaquette(Umu);
+      Real link_trace=WilsonLoops<Gimpl>::linkTrace(Umu); 
+      if( (orthog>=0) && (orthog<Nd) ){
+	std::cout << GridLogMessage << " Gauge fixing to Coulomb gauge time="<<orthog<< " plaq= "<<plaq<<" link trace = "<<link_trace<<  std::endl;
+      } else { 
+	std::cout << GridLogMessage << " Gauge fixing to Landau gauge plaq= "<<plaq<<" link trace = "<<link_trace<<  std::endl;
+      }
+    }
     for(int i=0;i<maxiter;i++){
 
       for(int mu=0;mu<Nd;mu++) U[mu]= PeekIndex<LorentzIndex>(Umu,mu);
 
       if ( Fourier==false ) { 
-	trG = SteepestDescentStep(U,xform,alpha,dmuAmu);
+	trG = SteepestDescentStep(U,xform,alpha,dmuAmu,orthog);
       } else { 
-	trG = FourierAccelSteepestDescentStep(U,xform,alpha,dmuAmu);
+	trG = FourierAccelSteepestDescentStep(U,xform,alpha,dmuAmu,orthog);
       }
+
+      //      std::cout << GridLogMessage << "trG   "<< trG<< std::endl;
+      //      std::cout << GridLogMessage << "xform "<< norm2(xform)<< std::endl;
+      //      std::cout << GridLogMessage << "dmuAmu "<< norm2(dmuAmu)<< std::endl;
+
       for(int mu=0;mu<Nd;mu++) PokeIndex<LorentzIndex>(Umu,U[mu],mu);
       // Monitor progress and convergence test 
       // infrequently to minimise cost overhead
@@ -106,14 +124,14 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
       }
     }
   };
-  static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu) {
+  static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) {
     GridBase *grid = U[0]._grid;
 
     std::vector<GaugeMat> A(Nd,grid);
     GaugeMat g(grid);
 
     GaugeLinkToLieAlgebraField(U,A);
-    ExpiAlphaDmuAmu(A,g,alpha,dmuAmu);
+    ExpiAlphaDmuAmu(A,g,alpha,dmuAmu,orthog);
 
 
     Real vol = grid->gSites();
@@ -125,7 +143,7 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
     return trG;
   }
 
-  static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu) {
+  static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) {
 
     GridBase *grid = U[0]._grid;
 
@@ -144,31 +162,41 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
 
     GaugeLinkToLieAlgebraField(U,A);
 
-    DmuAmu(A,dmuAmu);
+    DmuAmu(A,dmuAmu,orthog);
 
-    theFFT.FFT_all_dim(dmuAmu_p,dmuAmu,FFT::forward);
+    std::vector<int> mask(Nd,1);
+    for(int mu=0;mu<Nd;mu++) if (mu==orthog) mask[mu]=0;
+    theFFT.FFT_dim_mask(dmuAmu_p,dmuAmu,mask,FFT::forward);
 
     //////////////////////////////////
     // Work out Fp = psq_max/ psq...
+    // Avoid singularities in Fp
     //////////////////////////////////
     std::vector<int> latt_size = grid->GlobalDimensions();
     std::vector<int> coor(grid->_ndimension,0);
     for(int mu=0;mu<Nd;mu++) {
-
-      Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
-      LatticeCoordinate(pmu,mu);
-      pmu = TwoPiL * pmu ;
-      psq = psq + 4.0*sin(pmu*0.5)*sin(pmu*0.5); 
+      if ( mu != orthog ) { 
+	Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+	LatticeCoordinate(pmu,mu);
+	pmu = TwoPiL * pmu ;
+	psq = psq + 4.0*sin(pmu*0.5)*sin(pmu*0.5); 
+      }
     }
 
     Complex psqMax(16.0);
     Fp =  psqMax*one/psq;
 
-    pokeSite(TComplex(1.0),Fp,coor);
-
+    pokeSite(TComplex(16.0),Fp,coor);
+    if( (orthog>=0) && (orthog<Nd) ){
+      for(int t=0;t<grid->GlobalDimensions()[orthog];t++){
+	coor[orthog]=t;
+	pokeSite(TComplex(16.0),Fp,coor);
+      }
+    }
+    
     dmuAmu_p  = dmuAmu_p * Fp; 
 
-    theFFT.FFT_all_dim(dmuAmu,dmuAmu_p,FFT::backward);
+    theFFT.FFT_dim_mask(dmuAmu,dmuAmu_p,mask,FFT::backward);
 
     GaugeMat ciadmam(grid);
     Complex cialpha(0.0,-alpha);
@@ -183,11 +211,11 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
     return trG;
   }
 
-  static void ExpiAlphaDmuAmu(const std::vector<GaugeMat> &A,GaugeMat &g,Real & alpha, GaugeMat &dmuAmu) {
+  static void ExpiAlphaDmuAmu(const std::vector<GaugeMat> &A,GaugeMat &g,Real & alpha, GaugeMat &dmuAmu,int orthog) {
     GridBase *grid = g._grid;
     Complex cialpha(0.0,-alpha);
     GaugeMat ciadmam(grid);
-    DmuAmu(A,dmuAmu);
+    DmuAmu(A,dmuAmu,orthog);
     ciadmam = dmuAmu*cialpha;
     SU<Nc>::taExp(ciadmam,g);
   }  
diff --git a/tests/core/Test_fft_gfix.cc b/tests/core/Test_fft_gfix.cc
index 7fdf5f73..4eb6887d 100644
--- a/tests/core/Test_fft_gfix.cc
+++ b/tests/core/Test_fft_gfix.cc
@@ -60,6 +60,9 @@ int main (int argc, char ** argv)
   std::cout<< "* Testing we can gauge fix steep descent a RGT of Unit gauge    *" <<std::endl;
   std::cout<< "*****************************************************************" <<std::endl;
 
+  //  int coulomb_dir = -1;
+  int coulomb_dir = Nd-1;
+
   LatticeGaugeField   Umu(&GRID);
   LatticeGaugeField   Urnd(&GRID);
   LatticeGaugeField   Uorg(&GRID);
@@ -68,6 +71,7 @@ int main (int argc, char ** argv)
 
   LatticeColourMatrix   xform1(&GRID); // Gauge xform
   LatticeColourMatrix   xform2(&GRID); // Gauge xform
+  LatticeColourMatrix   xform3(&GRID); // Gauge xform
   
   SU3::ColdConfiguration(pRNG,Umu); // Unit gauge
   Uorg=Umu;
@@ -127,6 +131,22 @@ int main (int argc, char ** argv)
   plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
   std::cout << " Final plaquette "<<plaq << std::endl;
 
+  std::cout<< "*****************************************************************" <<std::endl;
+  std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge           *" <<std::endl;
+  std::cout<< "*****************************************************************" <<std::endl;
+
+  Umu=Urnd;
+  SU3::HotConfiguration(pRNG,Umu); // Unit gauge
+
+  plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
+  std::cout << " Initial plaquette "<<plaq << std::endl;
+
+  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
+
+  std::cout << Umu<<std::endl;
+
+  plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
+  std::cout << " Final plaquette "<<plaq << std::endl;
 
   Grid_finalize();
 }

From d14512ee03bf8bd6188f0ad26eac5cc59f19234e Mon Sep 17 00:00:00 2001
From: fionnoh <fionnoh@gmail.com>
Date: Fri, 17 May 2019 19:01:52 +0100
Subject: [PATCH 7/8] Exposed a coulomb/landau enum to the gauge fixing module

---
 Hadrons/Modules/MGauge/GaugeFix.hpp | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/Hadrons/Modules/MGauge/GaugeFix.hpp b/Hadrons/Modules/MGauge/GaugeFix.hpp
index 57b1e082..f5ebb39f 100644
--- a/Hadrons/Modules/MGauge/GaugeFix.hpp
+++ b/Hadrons/Modules/MGauge/GaugeFix.hpp
@@ -42,6 +42,8 @@ BEGIN_HADRONS_NAMESPACE
  ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MGauge)
 
+GRID_SERIALIZABLE_ENUM(Fix, undef, coulomb, Nd - 1, landau, -1);
+
 class GaugeFixPar: Serializable
 {
 public:
@@ -51,6 +53,7 @@ public:
                                     int, maxiter, 
                                     Real, Omega_tol, 
                                     Real, Phi_tol,
+                                    Fix,  gaugeFix,
                                     bool, Fourier);
 };
 
@@ -115,8 +118,8 @@ void TGaugeFix<GImpl>::execute(void)
 //Loads the gauge and fixes it
 {
     std::cout << "executing" << std::endl;
-    LOG(Message) << "Fixing the Gauge" << std::endl;
-    LOG(Message) << par().gauge << std::endl;
+    LOG(Message) << "Fixing the Gauge " << par().gauge << " using "
+                 << par().gaugeFix << " guage fixing. " << Nd - 1 << std::endl;
     auto &U     = envGet(GaugeField, par().gauge);
     auto &Umu   = envGet(GaugeField, getName());
     auto &xform = envGet(GaugeMat, getName()+"_xform");
@@ -126,9 +129,10 @@ void TGaugeFix<GImpl>::execute(void)
     int  maxiter   = par().maxiter;
     Real Omega_tol = par().Omega_tol;
     Real Phi_tol   = par().Phi_tol;
+    int  gaugeFix  = par().gaugeFix;
     bool Fourier   = par().Fourier;
     Umu = U;
-    FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier);
+    FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,gaugeFix);
     LOG(Message) << "Gauge Fixed" << std::endl;
 }
 

From dbd7f3f0fcbb0035c489e986d83fefd6b9a73b7c Mon Sep 17 00:00:00 2001
From: fionnoh <fionnoh@gmail.com>
Date: Fri, 17 May 2019 19:10:09 +0100
Subject: [PATCH 8/8] Added variables that were missing from wall source setup

---
 Hadrons/Modules/MSource/Wall.hpp | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/Hadrons/Modules/MSource/Wall.hpp b/Hadrons/Modules/MSource/Wall.hpp
index e43c61bb..945e2760 100644
--- a/Hadrons/Modules/MSource/Wall.hpp
+++ b/Hadrons/Modules/MSource/Wall.hpp
@@ -119,6 +119,9 @@ template <typename FImpl>
 void TWall<FImpl>::setup(void)
 {
     envCreateLat(PropagatorField, getName());
+    envCache(Lattice<iScalar<vInteger>>, tName_, 1, envGetGrid(LatticeComplex));
+    envCacheLat(LatticeComplex, momphName_);
+    envTmpLat(LatticeComplex, "coor");
 }
 
 // execution ///////////////////////////////////////////////////////////////////