Mres changes and gauge xform mat changes

Grid/qcd/action/fermion/FermionOperator.h

@@ -136,6 +136,7 @@ namespace Grid {
       //////////////////////////////////////////////////////////////////////
       // Conserved currents, either contract at sink or insert sequentially.
       //////////////////////////////////////////////////////////////////////
+
       virtual void ContractConservedCurrent(PropagatorField &q_in_1,
                                             PropagatorField &q_in_2,
                                             PropagatorField &q_out,
@@ -148,6 +149,12 @@ namespace Grid {
                                        unsigned int tmin, 
                                        unsigned int tmax,
                                        ComplexField &lattice_cmplx)=0;
+
+      // Only reimplemented in Wilson5D 
+      // Default to just a zero correlation function
+      virtual void ContractJ5q(FermionField &q_in   ,ComplexField &J5q) { J5q=zero; };
+      virtual void ContractJ5q(PropagatorField &q_in,ComplexField &J5q) { J5q=zero; };
+
       ///////////////////////////////////////////////
       // Physical field import/export
       ///////////////////////////////////////////////

Grid/qcd/action/fermion/WilsonFermion5D.cc

@@ -939,6 +939,75 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const Fe
     merge(qSiteRev, qSiteVec); \
 }
 
+//          psi = chiralProjectPlus(Result_s[Ls/2-1]);
+//          psi+= chiralProjectMinus(Result_s[Ls/2]);
+//         PJ5q+=localInnerProduct(psi,psi);
+
+template<class vobj> 
+Lattice<vobj> spProj5p(const Lattice<vobj> & in)
+{
+  GridBase *grid=in._grid;
+  Gamma G5(Gamma::Algebra::Gamma5);
+  Lattice<vobj> ret(grid);
+  parallel_for(int ss=0;ss<grid->oSites();ss++){
+    ret._odata[ss] = in._odata[ss] + G5*in._odata[ss];
+  }
+  return ret;
+}
+template<class vobj> 
+Lattice<vobj> spProj5m(const Lattice<vobj> & in)
+{
+  Gamma G5(Gamma::Algebra::Gamma5);
+  GridBase *grid=in._grid;
+  Lattice<vobj> ret(grid);
+  parallel_for(int ss=0;ss<grid->oSites();ss++){
+    ret._odata[ss] = in._odata[ss] - G5*in._odata[ss];
+  }
+  return ret;
+}
+
+template <class Impl>
+void WilsonFermion5D<Impl>::ContractJ5q(FermionField &q_in,ComplexField &J5q)
+{
+  conformable(GaugeGrid(), J5q._grid);
+  conformable(q_in._grid, FermionGrid());
+
+  // 4d field
+  int Ls = this->Ls;
+  FermionField psi(GaugeGrid());
+  FermionField p_plus (GaugeGrid());
+  FermionField p_minus(GaugeGrid());
+  FermionField p(GaugeGrid());
+
+  ExtractSlice(p_plus , q_in, Ls/2   , 0);
+  ExtractSlice(p_minus, q_in, Ls/2-1 , 0);
+  p_plus = spProj5p(p_plus );
+  p_minus= spProj5m(p_minus);
+  p=p_plus+p_minus;
+  J5q = localInnerProduct(p,p);
+}
+
+template <class Impl>
+void WilsonFermion5D<Impl>::ContractJ5q(PropagatorField &q_in,ComplexField &J5q)
+{
+  conformable(GaugeGrid(), J5q._grid);
+  conformable(q_in._grid, FermionGrid());
+
+  // 4d field
+  int Ls = this->Ls;
+  PropagatorField psi(GaugeGrid());
+  PropagatorField p_plus (GaugeGrid());
+  PropagatorField p_minus(GaugeGrid());
+  PropagatorField p(GaugeGrid());
+
+  ExtractSlice(p_plus , q_in, Ls/2   , 0);
+  ExtractSlice(p_minus, q_in, Ls/2-1 , 0);
+  p_plus = spProj5p(p_plus );
+  p_minus= spProj5m(p_minus);
+  p=p_plus+p_minus;
+  J5q = localInnerProduct(p,p);
+}
+
 template <class Impl>
 void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
                                                      PropagatorField &q_in_2,
@@ -949,6 +1018,7 @@ void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
     conformable(q_in_1._grid, FermionGrid());
     conformable(q_in_1._grid, q_in_2._grid);
     conformable(_FourDimGrid, q_out._grid);
+
     PropagatorField tmp1(FermionGrid()), tmp2(FermionGrid());
     unsigned int LLs = q_in_1._grid->_rdimensions[0];
     q_out = zero;
@@ -995,7 +1065,6 @@ void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
 }
 
 
-
 template <class Impl>
 void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
                                                 PropagatorField &q_out,

Grid/qcd/action/fermion/WilsonFermion5D.h

@@ -230,6 +230,10 @@ namespace QCD {
                              unsigned int tmin, 
                              unsigned int tmax,
 			     ComplexField &lattice_cmplx);
+
+    void ContractJ5q(PropagatorField &q_in,ComplexField &J5q);
+    void ContractJ5q(FermionField &q_in,ComplexField &J5q);
+
   };
 
 }}

Grid/qcd/utils/GaugeFix.h

@@ -53,21 +53,32 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
   }  
   static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false) {
     GridBase *grid = Umu._grid;
+    GaugeMat xform(grid);
+    SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier);
+  }
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false) {
+
+    GridBase *grid = Umu._grid;
 
     Real org_plaq      =WilsonLoops<Gimpl>::avgPlaquette(Umu);
     Real org_link_trace=WilsonLoops<Gimpl>::linkTrace(Umu); 
     Real old_trace = org_link_trace;
     Real trG;
+    
+    xform=1.0;
 
     std::vector<GaugeMat> U(Nd,grid);
-                 GaugeMat dmuAmu(grid);
+
+    GaugeMat dmuAmu(grid);
 
     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,alpha,dmuAmu);
+	trG = SteepestDescentStep(U,xform,alpha,dmuAmu);
       } else { 
-	trG = FourierAccelSteepestDescentStep(U,alpha,dmuAmu);
+	trG = FourierAccelSteepestDescentStep(U,xform,alpha,dmuAmu);
       }
       for(int mu=0;mu<Nd;mu++) PokeIndex<LorentzIndex>(Umu,U[mu],mu);
       // Monitor progress and convergence test 
@@ -84,7 +95,6 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
 	Real Phi  = 1.0 - old_trace / link_trace ;
 	Real Omega= 1.0 - trG;
 
-
 	std::cout << GridLogMessage << " Iteration "<<i<< " Phi= "<<Phi<< " Omega= " << Omega<< " trG " << trG <<std::endl;
 	if ( (Omega < Omega_tol) && ( ::fabs(Phi) < Phi_tol) ) {
 	  std::cout << GridLogMessage << "Converged ! "<<std::endl;
@@ -96,7 +106,7 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
       }
     }
   };
-  static Real SteepestDescentStep(std::vector<GaugeMat> &U,Real & alpha, GaugeMat & dmuAmu) {
+  static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu) {
     GridBase *grid = U[0]._grid;
 
     std::vector<GaugeMat> A(Nd,grid);
@@ -109,12 +119,13 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
     Real vol = grid->gSites();
     Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
 
+    xform = g*xform ;
     SU<Nc>::GaugeTransform(U,g);
 
     return trG;
   }
 
-  static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,Real & alpha, GaugeMat & dmuAmu) {
+  static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu) {
 
     GridBase *grid = U[0]._grid;
 
@@ -153,13 +164,6 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
     Complex psqMax(16.0);
     Fp =  psqMax*one/psq;
 
-    /*
-    static int once;
-    if ( once == 0 ) { 
-      std::cout << " Fp " << Fp <<std::endl;
-      once ++;
-      }*/
-
     pokeSite(TComplex(1.0),Fp,coor);
 
     dmuAmu_p  = dmuAmu_p * Fp; 
@@ -173,6 +177,7 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
 
     Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
 
+    xform = g*xform ;
     SU<Nc>::GaugeTransform(U,g);
 
     return trG;

Grid/qcd/utils/SUn.h

@@ -676,10 +676,18 @@ class SU {
     }
   }
 /*
- add GaugeTrans
+ * Fundamental rep gauge xform
-*/
+ */
-
+  template<typename Fundamental,typename GaugeMat>
-template<typename GaugeField,typename GaugeMat>
+  static void GaugeTransformFundamental( Fundamental &ferm, GaugeMat &g){
+    GridBase *grid = ferm._grid;
+    conformable(grid,g._grid);
+    ferm = g*ferm;
+  }
+/*
+ * Adjoint rep gauge xform
+ */
+  template<typename GaugeField,typename GaugeMat>
   static void GaugeTransform( GaugeField &Umu, GaugeMat &g){
     GridBase *grid = Umu._grid;
     conformable(grid,g._grid);

tests/core/Test_fft_gfix.cc

@@ -63,8 +63,11 @@ int main (int argc, char ** argv)
   LatticeGaugeField   Umu(&GRID);
   LatticeGaugeField   Urnd(&GRID);
   LatticeGaugeField   Uorg(&GRID);
+  LatticeGaugeField   Utmp(&GRID);
   LatticeColourMatrix   g(&GRID); // Gauge xform
 
+  LatticeColourMatrix   xform1(&GRID); // Gauge xform
+  LatticeColourMatrix   xform2(&GRID); // Gauge xform
   
   SU3::ColdConfiguration(pRNG,Umu); // Unit gauge
   Uorg=Umu;
@@ -78,7 +81,14 @@ int main (int argc, char ** argv)
   Real alpha=0.1;
 
   Umu = Urnd;
-  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,false);
+  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
+
+  // Check the gauge xform matrices
+  Utmp=Urnd;
+  SU<Nc>::GaugeTransform(Utmp,xform1);
+  Utmp = Utmp - Umu;
+  std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
+  
 
   plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
   std::cout << " Final plaquette "<<plaq << std::endl;
@@ -92,7 +102,13 @@ int main (int argc, char ** argv)
   std::cout<< "* Testing Fourier accelerated fixing                            *" <<std::endl;
   std::cout<< "*****************************************************************" <<std::endl;
   Umu=Urnd;
-  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
+  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
+
+  Utmp=Urnd;
+  SU<Nc>::GaugeTransform(Utmp,xform2);
+  Utmp = Utmp - Umu;
+  std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
+
 
   plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
   std::cout << " Final plaquette "<<plaq << std::endl;