From 0d9aa8722845c5b240445d93b009c2e9aa134dfd Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Thu, 22 Apr 2021 15:56:59 -0400
Subject: [PATCH] Reduce momentum to the GP plane

---
 tests/forces/Test_gpwilson_force.cc | 29 ++++++++++++++++++++++-------
 1 file changed, 22 insertions(+), 7 deletions(-)

diff --git a/tests/forces/Test_gpwilson_force.cc b/tests/forces/Test_gpwilson_force.cc
index d731f27a..65bd865e 100644
--- a/tests/forces/Test_gpwilson_force.cc
+++ b/tests/forces/Test_gpwilson_force.cc
@@ -64,7 +64,9 @@ int main (int argc, char ** argv)
   ////////////////////////////////////
   RealD mass=0.01; 
 
-  const int nu = 3;
+  const int nu = 1;
+  const int Lnu=latt_size[nu];
+
   std::vector<int> twists(Nd,0);  twists[nu] = 1;
   GparityWilsonFermionR::ImplParams params;  params.twists = twists;
   GparityWilsonFermionR Wil(U,*UGrid,*UrbGrid,mass,params);
@@ -84,20 +86,31 @@ int main (int argc, char ** argv)
   ////////////////////////////////////
   // Modify the gauge field a little 
   ////////////////////////////////////
-  RealD dt = 0.01;
+  RealD dt = 0.1;
 
   LatticeColourMatrix mommu(UGrid); 
+  LatticeColourMatrix zz(UGrid);
   LatticeColourMatrix forcemu(UGrid); 
   LatticeGaugeField mom(UGrid); 
   LatticeGaugeField Uprime(UGrid); 
 
+  
+  Lattice<iScalar<vInteger> > coor(UGrid);
+  LatticeCoordinate(coor,nu);
+  zz=Zero();
   for(int mu=0;mu<Nd;mu++){
 
     // Traceless antihermitian momentum; gaussian in lie alg
-    SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); 
-
+    if(mu==nu){
+      SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); 
+      mommu=where(coor==Lnu-1,mommu,zz);
+      //mommu=where(coor==0,mommu,zz);
+      //      mommu=where(coor==1,mommu,zz);
+    } else {
+      mommu=Zero();
+    }
     PokeIndex<LorentzIndex>(mom,mommu,mu);
-
+    
     // fourth order exponential approx
     autoView( mom_v, mom, CpuRead);
     autoView( U_v , U, CpuRead);
@@ -139,12 +152,14 @@ int main (int argc, char ** argv)
     // Update PF action density
     dS = dS+trace(mommu*forcemu)*dt;
   }
-
+  std::cout << "mommu"<<mommu<<std::endl;
+  std::cout << "dS" << dS<<std::endl;
+  
   ComplexD dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
-  std::cout << GridLogMessage << "dS      "<<Sprime-S<<std::endl;
+  std::cout << GridLogMessage << "Delta S "<<Sprime-S<<std::endl;
   std::cout << GridLogMessage << "predict dS    "<< dSpred <<std::endl;
 
   assert( fabs(real(Sprime-S-dSpred)) < 2.0 ) ;