From d7767a2a62508fae7450d69740f47cdaef915021 Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Fri, 10 Mar 2017 22:33:48 +0900
Subject: [PATCH] Few more tests

---
 lib/qcd/hmc/integrators/Integrator.h          | 21 ++++++++++---------
 .../hmc/integrators/Integrator_algorithm.h    |  3 ++-
 lib/qcd/utils/Metric.h                        | 10 ++++-----
 3 files changed, 18 insertions(+), 16 deletions(-)

diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index bba8b513..b162bcd2 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -125,10 +125,10 @@ class Integrator {
     Mom -= MomDer * ep;
 
     // Auxiliary fields
-    //P.update_auxiliary_momenta(ep*0.5);
-    //P.AuxiliaryFieldsDerivative(MomDer);
-    //Mom -= MomDer * ep;
-    //P.update_auxiliary_momenta(ep*0.5);
+    P.update_auxiliary_momenta(ep*0.5);
+    P.AuxiliaryFieldsDerivative(MomDer);
+    Mom -= MomDer * ep;
+    P.update_auxiliary_momenta(ep*0.5);
 
 
     for (int a = 0; a < as[level].actions.size(); ++a) {
@@ -193,9 +193,9 @@ class Integrator {
     }
 
     // Auxiliary fields
-    //P.update_auxiliary_momenta(ep*0.5);
-    //P.AuxiliaryFieldsDerivative(AuxDer);
-    //Msum += AuxDer;
+    P.update_auxiliary_momenta(ep*0.5);
+    P.AuxiliaryFieldsDerivative(AuxDer);
+    Msum += AuxDer;
     
 
     // Here run recursively
@@ -222,7 +222,7 @@ class Integrator {
     P.Mom = NewMom;
 
     // update the auxiliary fields momenta    
-    //P.update_auxiliary_momenta(ep*0.5);
+    P.update_auxiliary_momenta(ep*0.5);
   }
 
 
@@ -265,7 +265,7 @@ class Integrator {
     P.DerivativeP(Mom1); // first term in the derivative 
 
  
-    //P.update_auxiliary_fields(ep*0.5);
+    P.update_auxiliary_fields(ep*0.5);
 
 
     do {
@@ -293,7 +293,7 @@ class Integrator {
 
     U = NewU;
 
-    //P.update_auxiliary_fields(ep*0.5);
+    P.update_auxiliary_fields(ep*0.5);
 
   }
 
@@ -340,6 +340,7 @@ class Integrator {
 
   void reverse_momenta(){
     P.Mom *= -1.0;
+    P.AuxMom *= -1.0;
   }
 
   // to be used by the actionlevel class to iterate
diff --git a/lib/qcd/hmc/integrators/Integrator_algorithm.h b/lib/qcd/hmc/integrators/Integrator_algorithm.h
index 3e41a266..d61e28ee 100644
--- a/lib/qcd/hmc/integrators/Integrator_algorithm.h
+++ b/lib/qcd/hmc/integrators/Integrator_algorithm.h
@@ -344,7 +344,8 @@ class ImplicitLeapFrog : public Integrator<FieldImplementation, SmearingPolicy,
       if (last_step){
         this->update_P(U, level, eps / 2.0);
       } else {
-      this->implicit_update_P(U, level, eps, true);// intermediate step
+      this->implicit_update_P(U, level, eps, true);// works intermediate step
+      //  this->update_P(U, level, eps); // looks not reversible
       }
     }
   }
diff --git a/lib/qcd/utils/Metric.h b/lib/qcd/utils/Metric.h
index 8eeb5846..60a9bfc5 100644
--- a/lib/qcd/utils/Metric.h
+++ b/lib/qcd/utils/Metric.h
@@ -105,7 +105,7 @@ public:
     // Modify the distribution with the metric
     M.MSquareRoot(Mom);
 
-    if (0) {
+    if (1) {
       // Auxiliary momenta
       // do nothing if trivial, so hide in the metric
       MomentaField AuxMomTemp(Mom._grid);
@@ -132,7 +132,7 @@ public:
       Hloc += trace(Mom_mu * inv_mu);
     }
 
-    if (0) {
+    if (1) {
       // Auxiliary Fields
       // hide in the metric
       M.M(AuxMom, inv);
@@ -167,7 +167,7 @@ public:
 
   void AuxiliaryFieldsDerivative(MomentaField& der){
     der = zero;
-    if (0){
+    if (1){
     // Auxiliary fields
     MomentaField der_temp(der._grid);
     MomentaField X(der._grid);
@@ -194,13 +194,13 @@ public:
   }
 
   void update_auxiliary_momenta(RealD ep){
-    if(0){
+    if(1){
       AuxMom -= ep * AuxField;
     }
   }
 
   void update_auxiliary_fields(RealD ep){
-    if (0) {
+    if (1) {
       MomentaField tmp(AuxMom._grid);
       MomentaField tmp2(AuxMom._grid);
       M.M(AuxMom, tmp);