From f4ebea3381046026276864f3f908cb10b114d6a5 Mon Sep 17 00:00:00 2001 From: James Harrison Date: Mon, 14 Nov 2016 17:51:53 +0000 Subject: [PATCH] QedFVol: add functions for computing spatial and timelike Wilson loops --- programs/qed-fvol/WilsonLoops.h | 117 +++++++++++++++++++++++++++++--- programs/qed-fvol/qed-fvol.cc | 8 ++- 2 files changed, 114 insertions(+), 11 deletions(-) diff --git a/programs/qed-fvol/WilsonLoops.h b/programs/qed-fvol/WilsonLoops.h index c40fbaf3..98db6b7a 100644 --- a/programs/qed-fvol/WilsonLoops.h +++ b/programs/qed-fvol/WilsonLoops.h @@ -45,7 +45,7 @@ public: const std::vector &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 U(4, Umu._grid); - for (int mu = 0; mu < Nd; mu++) { + for (int mu = 0; mu < Umu._grid->_ndimension; mu++) { U[mu] = PeekIndex(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(); - double faces = (1.0 * Nd * (Nd - 1)) / 2.0; - return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME + Real vol = Umu._grid->gSites(); + Real faces = (1.0 * ndim * (ndim - 1)) / 2.0; + 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 &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 &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 U(4, Umu._grid); - for (int mu = 0; mu < Nd; mu++) { + for (int mu = 0; mu < Umu._grid->_ndimension; mu++) { U[mu] = PeekIndex(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 U(4, Umu._grid); + + for (int mu = 0; mu < Umu._grid->_ndimension; mu++) { + U[mu] = PeekIndex(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 U(4, Umu._grid); + + for (int mu = 0; mu < Umu._grid->_ndimension; mu++) { + U[mu] = PeekIndex(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(); - double faces = 1.0 * Nd * (Nd - 1); - return sumWl / vol / faces / Nc; // Nd , Nc dependent... FIXME + Real vol = Umu._grid->gSites(); + Real faces = 1.0 * ndim * (ndim - 1); + 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 } }; diff --git a/programs/qed-fvol/qed-fvol.cc b/programs/qed-fvol/qed-fvol.cc index 68705b8f..d026057e 100644 --- a/programs/qed-fvol/qed-fvol.cc +++ b/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::avgPlaquette(expA); avgWl2x2Aexp = NewWilsonLoops::avgWilsonLoop(expA, 2, 2); + avgWl2x2Aexp_time = NewWilsonLoops::avgTimelikeWilsonLoop(expA, 2, 2); + avgWl2x2Aexp_space = NewWilsonLoops::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