From f61c241c18e5e43628ab7b0a59d505f4fdd325ab Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Fri, 19 Jan 2018 09:52:45 +0000
Subject: [PATCH] Debugging shGordon, works

---
 lib/qcd/action/scalar/ScalarImpl.h     | 12 +++++++-----
 lib/qcd/action/scalar/shGordonAction.h |  8 ++++++--
 lib/qcd/hmc/integrators/Integrator.h   |  1 +
 3 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/lib/qcd/action/scalar/ScalarImpl.h b/lib/qcd/action/scalar/ScalarImpl.h
index f072960b..65acbb9e 100644
--- a/lib/qcd/action/scalar/ScalarImpl.h
+++ b/lib/qcd/action/scalar/ScalarImpl.h
@@ -29,9 +29,9 @@ class ScalarImplTypes {
     static inline Field projectForce(Field& P){return P;}
 
     static inline void update_field(Field& P, Field& U, double ep) {
-      std::cout << GridLogDebug << "P:\n" << P << std::endl;
+      //std::cout << GridLogDebug << "P:\n" << P << std::endl;
       U += P*ep;
-      std::cout << GridLogDebug << "U:\n" << U << std::endl;
+      //std::cout << GridLogDebug << "U:\n" << U << std::endl;
     }
 
     static inline RealD FieldSquareNorm(Field& U) {
@@ -39,15 +39,17 @@ class ScalarImplTypes {
     }
 
     static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
-      gaussian(pRNG, U);
+     random(pRNG, U);
     }
 
     static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
-      gaussian(pRNG, U);
+      random(pRNG, U);
+      U *= 0.01;
     }
 
     static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
-      U = 1.0;
+      U = 0.0;
+      //std::cout << GridLogDebug << "Initial U:\n" << U << std::endl;
     }
 
     static void MomentumSpacePropagator(Field &out, RealD m)
diff --git a/lib/qcd/action/scalar/shGordonAction.h b/lib/qcd/action/scalar/shGordonAction.h
index 64c94dd5..8eff60ab 100644
--- a/lib/qcd/action/scalar/shGordonAction.h
+++ b/lib/qcd/action/scalar/shGordonAction.h
@@ -60,12 +60,16 @@ class shGordonAction : public QCD::Action<typename Impl::Field> {
 
     virtual void deriv(const Field &phi,
                        Field &force) {
+        //std::cout << GridLogDebug << "Force total before :\n" << force << std::endl;
         Field tmp(phi._grid);
         tmp = 2.0*QCD::Nd*phi;
         for (int mu = 0; mu < QCD::Nd; mu++) tmp -= Cshift(phi, mu, 1) + Cshift(phi, mu, -1);
 
-        std::cout << GridLogDebug << "Phi norm : " << norm2(phi) << std::endl;
-        force+= tmp + 0.5*mass_square/g*(exp(g*phi) - exp(-g*phi));
+
+        //std::cout << GridLogDebug << "Phi norm : " << norm2(phi) << std::endl;
+        force += tmp + 0.5*mass_square/g*(exp(g*phi) - exp(-g*phi));
+        //std::cout << GridLogDebug << "Force tmp :\n" << tmp << std::endl;
+        //std::cout << GridLogDebug << "Force total after :\n" << force << std::endl;
     }
 };
 
diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index 9b99963d..a31d7e22 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -116,6 +116,7 @@ class Integrator {
 
    for (int a = 0; a < as[level].actions.size(); ++a) {
       Field force(U._grid);
+      force = zero;
       conformable(U._grid, Mom._grid);
       Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
       as[level].actions.at(a)->deriv(Us, force);  // deriv should NOT include Ta