From a80e43dbcf4b3106f75ee2c35ea56e6145fd4c9b Mon Sep 17 00:00:00 2001
From: James Harrison <jch1g10@soton.ac.uk>
Date: Wed, 11 Oct 2017 16:44:51 -0400
Subject: [PATCH] Added infinite-volume photon in Photon.h (not checked yet)

---
 lib/qcd/action/gauge/Photon.h | 112 +++++++++++++++++++++++++++++-----
 1 file changed, 96 insertions(+), 16 deletions(-)

diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h
index 7e21a1de..c8e3109a 100644
--- a/lib/qcd/action/gauge/Photon.h
+++ b/lib/qcd/action/gauge/Photon.h
@@ -58,7 +58,7 @@ namespace QCD{
   public:
     INHERIT_GIMPL_TYPES(Gimpl);
     GRID_SERIALIZABLE_ENUM(Gauge, undef, feynman, 1, coulomb, 2, landau, 3);
-    GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2);
+    GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2, qedInf, 3);
   public:
     Photon(Gauge gauge, ZmScheme zmScheme);
     virtual ~Photon(void) = default;
@@ -69,6 +69,7 @@ namespace QCD{
     void StochasticField(GaugeField &out, GridParallelRNG &rng,
                          const GaugeLinkField &weight);
   private:
+    void infVolPropagator(GaugeLinkField &out);
     void invKHatSquared(GaugeLinkField &out);
     void zmSub(GaugeLinkField &out);
   private:
@@ -95,6 +96,32 @@ namespace QCD{
     MomentumSpacePropagator(prop_k,in_k);
     theFFT.FFT_all_dim(out,prop_k,FFT::backward);
   }
+
+  template<class Gimpl>
+  void Photon<Gimpl>::infVolPropagator(GaugeLinkField &out)
+  {
+    GridBase           *grid = out._grid;
+    GaugeLinkField     xmu(grid), one(grid);
+    const unsigned int nd    = grid->_ndimension;
+    std::vector<int>   &l    = grid->_fdimensions;
+    std::vector<int>   x0(nd,0);
+    TComplex           Tone  = Complex(1.0,0.0);
+    TComplex           Tzero = Complex(0.0,0.0);
+    FFT                fft(grid);
+    
+    one = Complex(1.0,0.0);
+    out = zero;
+    for(int mu = 0; mu < nd; mu++)
+    {
+      LatticeCoordinate(xmu,mu);
+      xmu = where(xmu < latt_size[mu]/2, xmu, xmu-latt_size[mu]/2);
+      out = out + 4*M_PI*M_PI*xmu*xmu;
+    }
+    pokeSite(Tone, out, x0);
+    out = one/out;
+    pokeSite(Tzero, out, x0);
+    fft.FFT_all_dim(out, out, FFT::forward);
+  }
   
   template<class Gimpl>
   void Photon<Gimpl>::invKHatSquared(GaugeLinkField &out)
@@ -163,13 +190,28 @@ namespace QCD{
   void Photon<Gimpl>::MomentumSpacePropagator(const GaugeField &in,
                                                GaugeField &out)
   {
-    GridBase           *grid = out._grid;
-    LatticeComplex     k2Inv(grid);
+  GridBase           *grid = out._grid;
+    LatticeComplex     momProp(grid);
     
-    invKHatSquared(k2Inv);
-    zmSub(k2Inv);
+    switch (zmScheme_)
+    {
+      case ZmScheme::qedTL:
+      case ZmScheme::qedL:
+      {
+        invKHatSquared(momProp);
+        zmSub(momProp);
+        break;
+      }
+      case ZmScheme::qedInf:
+      {
+        infVolPropagator(momProp);
+        break;
+      }
+      default:
+        break;
+    }
     
-    out = in*k2Inv;
+    out = in*momProp;
   }
   
   template<class Gimpl>
@@ -179,14 +221,30 @@ namespace QCD{
     const unsigned int nd        = grid->_ndimension;
     std::vector<int>   latt_size = grid->_fdimensions;
     
-    Integer vol = 1;
-    for(int d = 0; d < nd; d++)
+    switch (zmScheme_)
     {
-      vol = vol * latt_size[d];
+      case ZmScheme::qedTL:
+      case ZmScheme::qedL:
+      {
+        Integer vol = 1;
+        for(int d = 0; d < nd; d++)
+        {
+          vol = vol * latt_size[d];
+        }
+        invKHatSquared(weight);
+        weight = sqrt(vol*real(weight));
+        zmSub(weight);
+        break;
+      }
+      case ZmScheme::qedInf:
+      {
+        infVolPropagator(weight);
+        weight = sqrt(real(weight));
+        break;
+      }
+      default:
+        break;
     }
-    invKHatSquared(weight);
-    weight = sqrt(vol*real(weight));
-    zmSub(weight);
   }
   
   template<class Gimpl>
@@ -209,12 +267,34 @@ namespace QCD{
     GaugeField         aTilde(grid);
     FFT                fft(grid);
     
-    for(int mu = 0; mu < nd; mu++)
+    switch (zmScheme_)
     {
-      gaussian(rng, r);
-      r = weight*r;
-      pokeLorentz(aTilde, r, mu);
+      case ZmScheme::qedTL:
+      case ZmScheme::qedL:
+      {
+        for(int mu = 0; mu < nd; mu++)
+        {
+          gaussian(rng, r);
+          r = weight*r;
+          pokeLorentz(aTilde, r, mu);
+        }
+        break;
+      }
+      case ZmScheme::qedInf:
+      {
+        Complex                    shift(1., 1.); // This needs to be a GaugeLink element?
+        for(int mu = 0; mu < nd; mu++)
+        {
+          bernoulli(rng, r);
+          r = weight*(2.*r - shift);
+          pokeLorentz(aTilde, r, mu);
+        }
+        break;
+      }
+      default:
+        break;
     }
+
     fft.FFT_all_dim(out, aTilde, FFT::backward);
     
     out = real(out);