From 048397d880f4c74c1cc18b7cfcb9565ff39dd6cc Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Sun, 24 Feb 2019 14:43:22 +0000
Subject: [PATCH] Default tau spacing should be longer c.f. Zbigniew Srocinsky
 thesis

---
 Grid/qcd/action/gauge/GaugeImplTypes.h        | 26 ++++++++++++++++---
 .../action/gauge/PlaqPlusRectangleAction.h    |  2 +-
 2 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h
index 9e3e0d68..818e5405 100644
--- a/Grid/qcd/action/gauge/GaugeImplTypes.h
+++ b/Grid/qcd/action/gauge/GaugeImplTypes.h
@@ -87,12 +87,32 @@ public:
   ///////////////////////////////////////////////////////////
   // Move these to another class
   // HMC auxiliary functions
-  static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) {
-    // specific for SU gauge fields
+  static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) 
+  {
+    // Zbigniew Srocinsky thesis:
+    //
+    // P(p) =  N \Prod_{x\mu}e^-{1/2 Tr (p^2_mux)}
+    // 
+    // p_x,mu = c_x,mu,a T_a
+    //
+    // Tr p^2 =  sum_a,x,mu 1/2 (c_x,mu,a)^2
+    //
+    // Which implies P(p) =  N \Prod_{x,\mu,a} e^-{1/4 c_xmua^2  }
+    //
+    //                    =  N \Prod_{x,\mu,a} e^-{1/2 (c_xmua/sqrt{2})^2  }
+    // 
+    // Expect c' = cxmua/sqrt(2) to be a unit variance gaussian.
+    //
+    // Expect cxmua variance sqrt(2).
+    //
+    // Must scale the momentum by sqrt(2) up to invoke CPS and UKQCD conventions
+    //
     LinkField Pmu(P._grid);
-    Pmu = zero;
+    Pmu = Zero();
     for (int mu = 0; mu < Nd; mu++) {
       SU<Nrepresentation>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
+      RealD scale = ::sqrt(2) ;
+      Pmu = Pmu*scale;
       PokeIndex<LorentzIndex>(P, Pmu, mu);
     }
   }
diff --git a/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h b/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h
index 5bfd39b2..bdbc8479 100644
--- a/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h
+++ b/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h
@@ -75,7 +75,7 @@ namespace Grid{
       virtual void deriv(const GaugeField &Umu,GaugeField & dSdU) {
 	//extend Ta to include Lorentz indexes
 	RealD factor_p = c_plaq/RealD(Nc)*0.5;
-	RealD factor_r =   c_rect/RealD(Nc)*0.5;
+	RealD factor_r = c_rect/RealD(Nc)*0.5;
 
 	GridBase *grid = Umu._grid;