From 391fd9cc6a4c4d826eb71c2d8d761afaea3b7a69 Mon Sep 17 00:00:00 2001
From: david clarke <clarke.davida@gmail.com>
Date: Tue, 17 Oct 2023 14:57:15 -0600
Subject: [PATCH] try lepage term

---
 Grid/qcd/smearing/HISQSmearing.h | 52 +++++++++++++++++++++++++-------
 examples/Example_plaquette.cc    | 14 +++++----
 tests/smearing/Test_fatLinks.cc  |  4 +--
 3 files changed, 51 insertions(+), 19 deletions(-)

diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h
index 432184e0..44e14e85 100644
--- a/Grid/qcd/smearing/HISQSmearing.h
+++ b/Grid/qcd/smearing/HISQSmearing.h
@@ -81,8 +81,8 @@ struct SmearingParameters{
 
 
 /*!  @brief create fat links from link variables */
-template<class LGF> 
-class Smear_HISQ_fat {
+template<class LGF, class Gimpl> 
+class Smear_HISQ_fat : public Gimpl {
 
 private:
     GridCartesian* const _grid;
@@ -90,6 +90,8 @@ private:
 
 public:
 
+    INHERIT_GIMPL_TYPES(Gimpl);
+
     // Don't allow default values here.
     Smear_HISQ_fat(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) 
         : _grid(grid), 
@@ -141,6 +143,7 @@ public:
         // This is where contributions from the smearing get added together
         Ughost_fat=Zero();
 
+        // This loop handles 3-, 5-, and 7-link constructs, minus Lepage and Naik.
         for(int mu=0;mu<Nd;mu++) {
 
             // Create the accessors
@@ -161,8 +164,7 @@ public:
             stencilElement SE0, SE1, SE2, SE3, SE4;
             U3matrix U0, U1, U2, U3, U4, U5, W;
 
-            // 3-link
-            for(int site=0;site<U_v.size();site++){
+            for(int site=0;site<U_v.size();site++){ // ----------- 3-link
                 for(int nu=0;nu<Nd;nu++) {
                     if(nu==mu) continue;
                     int s = stencilIndex(mu,nu);
@@ -195,8 +197,7 @@ public:
                 }
             }
 
-            // 5-link
-            for(int site=0;site<U_v.size();site++){
+            for(int site=0;site<U_v.size();site++){ // ----------- 5-link
                 int sigmaIndex = 0;
                 for(int nu=0;nu<Nd;nu++) {
                     if(nu==mu) continue;
@@ -210,8 +211,6 @@ public:
                         SE3 = gStencil.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
                         SE4 = gStencil.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
 
-                        // gpermutes will be replaced with single line of code, combines load and permute
-                        // into one step. still in pull request stage 
                         U0 = coalescedReadGeneralPermute(      U_v[x_p_mu     ](nu ),SE0->_permute,Nd);
                         U1 = coalescedReadGeneralPermute(U_3link_v[x_p_nu     ](rho),SE1->_permute,Nd);
                         U2 = coalescedReadGeneralPermute(      U_v[x          ](nu ),SE2->_permute,Nd);
@@ -234,8 +233,7 @@ public:
                 }
             }
 
-            // 7-link
-            for(int site=0;site<U_v.size();site++){
+            for(int site=0;site<U_v.size();site++){ // ----------- 7-link
                 int sigmaIndex = 0;
                 for(int nu=0;nu<Nd;nu++) {
                     if(nu==mu) continue;
@@ -275,9 +273,41 @@ public:
 
         } // end mu loop
 
+        // c1, c3, c5, c7 construct contributions
         u_smr = Ghost.Extract(Ughost_fat) + lt.c_1*u_thin;
-    };
 
+        // Load up U and V std::vectors to access thin and smeared links.
+        std::vector<LatticeColourMatrix> U(Nd, u_thin.Grid());
+        std::vector<LatticeColourMatrix> V(Nd, u_smr.Grid());
+        for (int mu = 0; mu < Nd; mu++) {
+            U[mu] = PeekIndex<LorentzIndex>(u_thin, mu);
+            V[mu] = PeekIndex<LorentzIndex>(u_smr, mu);
+        }
+
+        // Compute LePage term from U_thin:
+        for(int mu=0;mu<Nd;mu++) {
+            for (int nu_h=1;nu_h<Nd;nu_h++) {
+                int nu=(mu+nu_h)%Nd;
+                                // nu, nu, mu, Back(nu), Back(nu)
+                V[mu] = V[mu] + lt.c_lp*Gimpl::CovShiftForward(U[nu],nu,
+                                          Gimpl::CovShiftForward(U[nu],nu,
+                                            Gimpl::CovShiftForward(U[mu],mu,
+                                              Gimpl::CovShiftBackward(U[nu],nu,
+                                                Gimpl::CovShiftIdentityBackward(U[nu],nu)))))
+                                // Back(nu), Back(nu), mu, nu, nu
+                              + lt.c_lp*Gimpl::CovShiftBackward(U[nu],nu,
+                                          Gimpl::CovShiftBackward(U[nu],nu,
+                                            Gimpl::CovShiftForward(U[mu],mu,
+                                              Gimpl::CovShiftForward(U[nu],nu,
+                                                Gimpl::CovShiftIdentityForward(U[nu],nu)))));
+            }
+        }
+
+        // Put V back into u_smr.
+        for (int mu = 0; mu < Nd; mu++) {
+            PokeIndex<LorentzIndex>(u_smr, V[mu], mu);
+        }
+    };
 
 //    void derivative(const GaugeField& Gauge) const {
 //    };
diff --git a/examples/Example_plaquette.cc b/examples/Example_plaquette.cc
index 17de4762..faf17d82 100644
--- a/examples/Example_plaquette.cc
+++ b/examples/Example_plaquette.cc
@@ -29,7 +29,6 @@ public:
     // Gimpl seems to be an arbitrary class. Within this class, it is expected that certain types are
     // already defined, things like Scalar and Field. This macro includes a bunch of #typedefs that
     // implement this equivalence at compile time.
-    // WARNING: The first time you include this or take it out, the compile time will increase a lot.
     INHERIT_GIMPL_TYPES(Gimpl);
 
     // Some example Gimpls can be found in GaugeImplementations.h, at the bottom. These are in turn built
@@ -53,12 +52,15 @@ public:
 
     // U_mu_nu(x)
     static void dirPlaquette(GaugeMat &plaq, const std::vector<GaugeMat> &U, const int mu, const int nu) {
-        // These CovShift calls seem to carry out the multiplication already. A positive shift moves the lattice 
-        // site x_mu = 1 in the RHS to x_mu = 0 in the result.
+        // Calls like CovShiftForward and CovShiftBackward have 3 arguments, and they multiply together
+        // the first and last argument. (Second arg gives the shift direction.) The CovShiftIdentityBackward
+        // has meanwhile only two arguments; it just returns the shifted (adjoint since backward) link. 
         plaq = Gimpl::CovShiftForward(U[mu],mu,
-                    Gimpl::CovShiftForward(U[nu],nu,
-                        Gimpl::CovShiftBackward(U[mu],mu,
-                            Gimpl::CovShiftIdentityBackward(U[nu], nu))));
+                   // Means Link*Cshift(field,mu,1), arguments are Link, mu, field in that order.
+                   Gimpl::CovShiftForward(U[nu],nu,
+                       Gimpl::CovShiftBackward(U[mu],mu,
+                           // This means Cshift(adj(Link), mu, -1)
+                           Gimpl::CovShiftIdentityBackward(U[nu], nu))));
     }
 
     // tr U_mu_nu(x)
diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc
index f5c7b5ca..1b24c7ca 100644
--- a/tests/smearing/Test_fatLinks.cc
+++ b/tests/smearing/Test_fatLinks.cc
@@ -90,14 +90,14 @@ int main (int argc, char** argv) {
     NerscIO::readConfiguration(Umu, header, conf_in);
 
     // Smear Umu and store result in U_smr
-    Smear_HISQ_fat<LatticeGaugeField> hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,0.);
+    Smear_HISQ_fat<LatticeGaugeField,PeriodicGimplD> hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,-1/8.);
     hisq_fat.smear(U_smr,Umu);
 
     NerscIO::writeConfiguration(U_smr,conf_out,"HISQ");
 
     // Test a C-style instantiation 
     double path_coeff[6] = {1, 2, 3, 4, 5, 6};
-    Smear_HISQ_fat<LatticeGaugeField> hisq_fat_Cstyle(&GRID,path_coeff);
+    Smear_HISQ_fat<LatticeGaugeField,PeriodicGimplD> hisq_fat_Cstyle(&GRID,path_coeff);
 
     // Make sure result doesn't change w.r.t. a trusted lattice
     NerscIO::readConfiguration(Umu, header, "nersc.l8t4b3360.357link.control");