QedFVol: add functions for computing spatial and timelike Wilson loops

programs/qed-fvol/WilsonLoops.h

@@ -45,7 +45,7 @@ public:
                             const std::vector<GaugeMat> &U) {
     LatticeComplex sitePlaq(U[0]._grid);
     Plaq = zero;
-    for (int mu = 1; mu < Nd; mu++) {
+    for (int mu = 1; mu < U[0]._grid->_ndimension; mu++) {
       for (int nu = 0; nu < mu; nu++) {
         traceDirPlaquette(sitePlaq, U, mu, nu);
         Plaq = Plaq + sitePlaq;
@@ -58,7 +58,7 @@ public:
   static Real sumPlaquette(const GaugeLorentz &Umu) {
     std::vector<GaugeMat> U(4, Umu._grid);
 
-    for (int mu = 0; mu < Nd; mu++) {
+    for (int mu = 0; mu < Umu._grid->_ndimension; mu++) {
       U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
     }
 
@@ -74,10 +74,11 @@ public:
   // average over all x,y,z,t and over all planes of plaquette
   //////////////////////////////////////////////////
   static Real avgPlaquette(const GaugeLorentz &Umu) {
+    int ndim = Umu._grid->_ndimension;
     Real sumplaq = sumPlaquette(Umu);
-    double vol = Umu._grid->gSites();
+    Real vol = Umu._grid->gSites();
-    double faces = (1.0 * Nd * (Nd - 1)) / 2.0;
+    Real faces = (1.0 * ndim * (ndim - 1)) / 2.0;
-    return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME
+    return sumplaq / vol / faces / Nc; // Nc dependent... FIXME
   }
 
   //////////////////////////////////////////////////
@@ -123,7 +124,42 @@ public:
                             const int R1, const int R2) {
     LatticeComplex siteWl(U[0]._grid);
     Wl = zero;
-    for (int mu = 1; mu < Nd; mu++) {
+    for (int mu = 1; mu < U[0]._grid->_ndimension; mu++) {
+      for (int nu = 0; nu < mu; nu++) {
+        traceWilsonLoop(siteWl, U, R1, R2, mu, nu);
+        Wl = Wl + siteWl;
+        traceWilsonLoop(siteWl, U, R2, R1, mu, nu);
+        Wl = Wl + siteWl;
+      }
+    }
+  }
+  //////////////////////////////////////////////////
+  // sum over planes of Wilson loop with length R1
+  // in the time direction
+  //////////////////////////////////////////////////
+  static void siteTimelikeWilsonLoop(LatticeComplex &Wl,
+                            const std::vector<GaugeMat> &U,
+                            const int R1, const int R2) {
+    LatticeComplex siteWl(U[0]._grid);
+
+    int ndim = U[0]._grid->_ndimension;
+
+    Wl = zero;
+    for (int nu = 0; nu < ndim - 1; nu++) {
+      traceWilsonLoop(siteWl, U, R1, R2, ndim-1, nu);
+      Wl = Wl + siteWl;
+    }
+  }
+  //////////////////////////////////////////////////
+  // sum Wilson loop over all planes orthogonal to the time direction
+  //////////////////////////////////////////////////
+  static void siteSpatialWilsonLoop(LatticeComplex &Wl,
+                            const std::vector<GaugeMat> &U,
+                            const int R1, const int R2) {
+    LatticeComplex siteWl(U[0]._grid);
+
+    Wl = zero;
+    for (int mu = 1; mu < U[0]._grid->_ndimension - 1; mu++) {
       for (int nu = 0; nu < mu; nu++) {
         traceWilsonLoop(siteWl, U, R1, R2, mu, nu);
         Wl = Wl + siteWl;
@@ -139,7 +175,7 @@ public:
                             const int R1, const int R2) {
     std::vector<GaugeMat> U(4, Umu._grid);
 
-    for (int mu = 0; mu < Nd; mu++) {
+    for (int mu = 0; mu < Umu._grid->_ndimension; mu++) {
       U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
     }
 
@@ -152,14 +188,75 @@ public:
     return p.real();
   }
   //////////////////////////////////////////////////
+  // sum over all x,y,z,t and over all planes of timelike Wilson loop
+  //////////////////////////////////////////////////
+  static Real sumTimelikeWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    std::vector<GaugeMat> U(4, Umu._grid);
+
+    for (int mu = 0; mu < Umu._grid->_ndimension; mu++) {
+      U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+    }
+
+    LatticeComplex Wl(Umu._grid);
+
+    siteTimelikeWilsonLoop(Wl, U, R1, R2);
+
+    TComplex Tp = sum(Wl);
+    Complex p = TensorRemove(Tp);
+    return p.real();
+  }
+  //////////////////////////////////////////////////
+  // sum over all x,y,z,t and over all planes of spatial Wilson loop
+  //////////////////////////////////////////////////
+  static Real sumSpatialWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    std::vector<GaugeMat> U(4, Umu._grid);
+
+    for (int mu = 0; mu < Umu._grid->_ndimension; mu++) {
+      U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+    }
+
+    LatticeComplex Wl(Umu._grid);
+
+    siteSpatialWilsonLoop(Wl, U, R1, R2);
+
+    TComplex Tp = sum(Wl);
+    Complex p = TensorRemove(Tp);
+    return p.real();
+  }
+  //////////////////////////////////////////////////
   // average over all x,y,z,t and over all planes of Wilson loop
   //////////////////////////////////////////////////
   static Real avgWilsonLoop(const GaugeLorentz &Umu,
                             const int R1, const int R2) {
+    int ndim = Umu._grid->_ndimension;
     Real sumWl = sumWilsonLoop(Umu, R1, R2);
-    double vol = Umu._grid->gSites();
+    Real vol = Umu._grid->gSites();
-    double faces = 1.0 * Nd * (Nd - 1);
+    Real faces = 1.0 * ndim * (ndim - 1);
-    return sumWl / vol / faces / Nc; // Nd , Nc dependent... FIXME
+    return sumWl / vol / faces / Nc; // Nc dependent... FIXME
+  }
+  //////////////////////////////////////////////////
+  // average over all x,y,z,t and over all planes of timelike Wilson loop
+  //////////////////////////////////////////////////
+  static Real avgTimelikeWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    int ndim = Umu._grid->_ndimension;
+    Real sumWl = sumTimelikeWilsonLoop(Umu, R1, R2);
+    Real vol = Umu._grid->gSites();
+    Real faces = 1.0 * (ndim - 1);
+    return sumWl / vol / faces / Nc; // Nc dependent... FIXME
+  }
+  //////////////////////////////////////////////////
+  // average over all x,y,z,t and over all planes of spatial Wilson loop
+  //////////////////////////////////////////////////
+  static Real avgSpatialWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    int ndim = Umu._grid->_ndimension;
+    Real sumWl = sumSpatialWilsonLoop(Umu, R1, R2);
+    Real vol = Umu._grid->gSites();
+    Real faces = 1.0 * (ndim - 1) * (ndim - 2);
+    return sumWl / vol / faces / Nc; // Nc dependent... FIXME
   }
 };
 

programs/qed-fvol/qed-fvol.cc

@@ -64,7 +64,7 @@ int main(int argc, char *argv[])
     EmField          a(&grid);
     EmField          expA(&grid);
 
-    Real avgPlaqAexp, avgWl2x2Aexp;
+    Real avgPlaqAexp, avgWl2x2Aexp, avgWl2x2Aexp_time, avgWl2x2Aexp_space;
 
     pRNG.SeedRandomDevice();
     photon.StochasticField(a, pRNG);
@@ -117,14 +117,20 @@ int main(int argc, char *argv[])
     // Calculate plaquette from exponentiated photon field
     avgPlaqAexp = NewWilsonLoops<QedPeriodicGimplR>::avgPlaquette(expA);
     avgWl2x2Aexp = NewWilsonLoops<QedPeriodicGimplR>::avgWilsonLoop(expA, 2, 2);
+    avgWl2x2Aexp_time = NewWilsonLoops<QedPeriodicGimplR>::avgTimelikeWilsonLoop(expA, 2, 2);
+    avgWl2x2Aexp_space = NewWilsonLoops<QedPeriodicGimplR>::avgSpatialWilsonLoop(expA, 2, 2);
 
     avgPlaqAexp = avgPlaqAexp*3;
     avgWl2x2Aexp = avgWl2x2Aexp*3;
+    avgWl2x2Aexp_time = avgWl2x2Aexp_time*3;
+    avgWl2x2Aexp_space = avgWl2x2Aexp_space*3;
 
     LOG(Message) << "Plaquette average (from A): " << avgPlaqA << std::endl;
     LOG(Message) << "Plaquette average (from exp(A)): " << avgPlaqAexp << std::endl;
     LOG(Message) << "2x2 Wilson Loop average (from A): " << avgWlA << std::endl;
     LOG(Message) << "2x2 Wilson Loop average (from exp(A)): " << avgWl2x2Aexp << std::endl;
+    LOG(Message) << "2x2 Wilson Loop timelike average (from exp(A)): " << avgWl2x2Aexp_time << std::endl;
+    LOG(Message) << "2x2 Wilson Loop spatial average (from exp(A)): " << avgWl2x2Aexp_space << std::endl;
     LOG(Message) << "Plaquette (one site): " << plaqsite / faces << std::endl;
 
     // epilogue