From cf167d0cd1c561bed3557eaf89350b8d8eb8d9b1 Mon Sep 17 00:00:00 2001
From: James Harrison <jch1g10@soton.ac.uk>
Date: Mon, 14 Nov 2016 17:02:29 +0000
Subject: [PATCH] QedFVol: implement exponentiation of photon field

---
 programs/qed-fvol/WilsonLoops.h | 19 ++++++++++---------
 programs/qed-fvol/qed-fvol.cc   | 32 +++++++++++++++++++++++++-------
 2 files changed, 35 insertions(+), 16 deletions(-)

diff --git a/programs/qed-fvol/WilsonLoops.h b/programs/qed-fvol/WilsonLoops.h
index 610fdc7b..c40fbaf3 100644
--- a/programs/qed-fvol/WilsonLoops.h
+++ b/programs/qed-fvol/WilsonLoops.h
@@ -5,7 +5,7 @@
 
 BEGIN_QEDFVOL_NAMESPACE
 
-template <class Gimpl> class WilsonLoops : public Gimpl {
+template <class Gimpl> class NewWilsonLoops : public Gimpl {
 public:
   INHERIT_GIMPL_TYPES(Gimpl);
 
@@ -55,7 +55,7 @@ public:
   //////////////////////////////////////////////////
   // sum over all x,y,z,t and over all planes of plaquette
   //////////////////////////////////////////////////
-  static RealD sumPlaquette(const GaugeLorentz &Umu) {
+  static Real sumPlaquette(const GaugeLorentz &Umu) {
     std::vector<GaugeMat> U(4, Umu._grid);
 
     for (int mu = 0; mu < Nd; mu++) {
@@ -73,8 +73,8 @@ public:
   //////////////////////////////////////////////////
   // average over all x,y,z,t and over all planes of plaquette
   //////////////////////////////////////////////////
-  static RealD avgPlaquette(const GaugeLorentz &Umu) {
-    RealD sumplaq = sumPlaquette(Umu);
+  static Real avgPlaquette(const GaugeLorentz &Umu) {
+    Real sumplaq = sumPlaquette(Umu);
     double vol = Umu._grid->gSites();
     double faces = (1.0 * Nd * (Nd - 1)) / 2.0;
     return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME
@@ -112,14 +112,14 @@ public:
                                 const int Rmu, const int Rnu,
                                 const int mu, const int nu) {
     GaugeMat sp(U[0]._grid);
-    WilsonLoop(sp, U, Rmu, Rnu, mu, nu);
+    wilsonLoop(sp, U, Rmu, Rnu, mu, nu);
     wl = trace(sp);
   }
   //////////////////////////////////////////////////
   // sum over all planes of Wilson loop
   //////////////////////////////////////////////////
   static void siteWilsonLoop(LatticeComplex &Wl,
-                            const std::vector<GaugeMat> &U
+                            const std::vector<GaugeMat> &U,
                             const int R1, const int R2) {
     LatticeComplex siteWl(U[0]._grid);
     Wl = zero;
@@ -135,7 +135,7 @@ public:
   //////////////////////////////////////////////////
   // sum over all x,y,z,t and over all planes of Wilson loop
   //////////////////////////////////////////////////
-  static RealD sumWilsonLoop(const GaugeLorentz &Umu,
+  static Real sumWilsonLoop(const GaugeLorentz &Umu,
                             const int R1, const int R2) {
     std::vector<GaugeMat> U(4, Umu._grid);
 
@@ -154,13 +154,14 @@ public:
   //////////////////////////////////////////////////
   // average over all x,y,z,t and over all planes of Wilson loop
   //////////////////////////////////////////////////
-  static RealD avgPlaquette(const GaugeLorentz &Umu,
+  static Real avgWilsonLoop(const GaugeLorentz &Umu,
                             const int R1, const int R2) {
-    RealD sumWl = sumWilsonLoop(Umu);
+    Real sumWl = sumWilsonLoop(Umu, R1, R2);
     double vol = Umu._grid->gSites();
     double faces = 1.0 * Nd * (Nd - 1);
     return sumWl / vol / faces / Nc; // Nd , Nc dependent... FIXME
   }
+};
 
 END_QEDFVOL_NAMESPACE
 
diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc
index fd780edf..68705b8f 100644
--- a/programs/qed-fvol/qed-fvol.cc
+++ b/programs/qed-fvol/qed-fvol.cc
@@ -1,4 +1,5 @@
 #include <Global.hpp>
+#include <WilsonLoops.h>
 
 using namespace Grid;
 using namespace QCD;
@@ -24,10 +25,11 @@ public:
   typedef Lattice<SiteGaugeField> GaugeField;
 };
 
-typedef QedGimpl<vComplex>      QedGimplR;
-typedef Photon<QedGimplR>       PhotonR;
-typedef PhotonR::GaugeField     EmField;
-typedef PhotonR::GaugeLinkField EmComp;
+typedef QedGimpl<vComplex>              QedGimplR;
+typedef PeriodicGaugeImpl<QedGimplR>    QedPeriodicGimplR;
+typedef Photon<QedGimplR>               PhotonR;
+typedef PhotonR::GaugeField             EmField;
+typedef PhotonR::GaugeLinkField         EmComp;
 
 int main(int argc, char *argv[])
 {
@@ -60,11 +62,18 @@ int main(int argc, char *argv[])
     PhotonR          photon(PhotonR::Gauge::Feynman,
                             PhotonR::ZmScheme::QedL);
     EmField          a(&grid);
+    EmField          expA(&grid);
+
+    Real avgPlaqAexp, avgWl2x2Aexp;
 
     pRNG.SeedRandomDevice();
     photon.StochasticField(a, pRNG);
 
-    // Calculate log of plaquette
+    // Exponentiate photon field
+    Complex imag_unit(0, 1);
+    expA = exp(imag_unit*0.5*(a+conjugate(a)));
+
+    // Calculate plaquette from photon field
     EmComp              plaqA(&grid);
     EmComp              wlA(&grid);
     EmComp              tmp(&grid);
@@ -105,8 +114,17 @@ int main(int argc, char *argv[])
     peekSite(tplaqsite, plaqtrace, site0);
     Complex plaqsite = TensorRemove(tplaqsite);
 
-    LOG(Message) << "Plaquette average: " << avgPlaqA << std::endl;
-    LOG(Message) << "2x2 Wilson Loop average: " << avgWlA << std::endl;
+    // Calculate plaquette from exponentiated photon field
+    avgPlaqAexp = NewWilsonLoops<QedPeriodicGimplR>::avgPlaquette(expA);
+    avgWl2x2Aexp = NewWilsonLoops<QedPeriodicGimplR>::avgWilsonLoop(expA, 2, 2);
+
+    avgPlaqAexp = avgPlaqAexp*3;
+    avgWl2x2Aexp = avgWl2x2Aexp*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) << "Plaquette (one site): " << plaqsite / faces << std::endl;
 
     // epilogue