From ec035983fd423780af71e1e8ef0107a60715b4ea Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Wed, 1 Mar 2017 11:56:35 +0000
Subject: [PATCH] Fixing the implicit integration

---
 lib/qcd/hmc/GenericHMCrunner.h             |  2 +-
 lib/qcd/hmc/integrators/Integrator.h       | 41 +++++++++++++---------
 tests/hmc/Test_hmc_WilsonGauge_Implicit.cc |  8 ++---
 3 files changed, 30 insertions(+), 21 deletions(-)

diff --git a/lib/qcd/hmc/GenericHMCrunner.h b/lib/qcd/hmc/GenericHMCrunner.h
index f8a27f69..d4626ed6 100644
--- a/lib/qcd/hmc/GenericHMCrunner.h
+++ b/lib/qcd/hmc/GenericHMCrunner.h
@@ -142,7 +142,7 @@ class HMCWrapperTemplate: public HMCRunnerBase<ReaderClass> {
     //TrivialMetric<typename Implementation::Field> Mtr;
     ConjugateGradient<LatticeGaugeField> CG(1.0e-8,10000);
     LaplacianParams LapPar(0.0001, 1.0, 1000, 1e-8, 12, 64);
-    RealD Kappa = 0.9;
+    RealD Kappa = 0.6;
 
     // Better to pass the generalised momenta to the integrator
     LaplacianAdjointField<PeriodicGimplR> Laplacian(UGrid, CG, LapPar, Kappa);
diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index 15861215..bba8b513 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -114,6 +114,23 @@ class Integrator {
   void update_P(MomentaField& Mom, Field& U, int level, double ep) {
     // input U actually not used in the fundamental case
     // Fundamental updates, include smearing
+
+    // Generalised momenta  
+    // Derivative of the kinetic term must be computed before
+    // Mom is the momenta and gets updated by the 
+    // actions derivatives
+    MomentaField MomDer(P.Mom._grid);
+    P.M.ImportGauge(U);
+    P.DerivativeU(P.Mom, MomDer);
+    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);
+
+
     for (int a = 0; a < as[level].actions.size(); ++a) {
       Field force(U._grid);
       conformable(U._grid, Mom._grid);
@@ -128,17 +145,7 @@ class Integrator {
       Mom -= force * ep; 
     }
 
-    // Generalised momenta
-    MomentaField MomDer(P.Mom._grid);
-    P.M.ImportGauge(U);
-    P.DerivativeU(P.Mom, MomDer);
-    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);
 
     // Force from the other representations
     as[level].apply(update_P_hireps, Representations, Mom, U, ep);
@@ -179,9 +186,11 @@ class Integrator {
     MomentaField AuxDer(P.Mom._grid);
     MomDer1 = zero;
     MomentaField diff(P.Mom._grid);
-
-    if (intermediate)
+    double factor = 2.0;
+    if (intermediate){
       P.DerivativeU(P.Mom, MomDer1);
+      factor = 1.0;
+    }
 
     // Auxiliary fields
     //P.update_auxiliary_momenta(ep*0.5);
@@ -202,7 +211,7 @@ class Integrator {
       std::cout << GridLogIntegrator << "|Force| laplacian site average: " << force_abs
                 << std::endl;
 
-      NewMom = P.Mom - ep* 0.5 * (2.0*Msum + MomDer + MomDer1);// simplify
+      NewMom = P.Mom - ep* 0.5 * (2.0*Msum + factor*MomDer + MomDer1);// simplify
       diff = NewMom - OldMom;
       counter++;
       RelativeError = std::sqrt(norm2(diff))/std::sqrt(norm2(NewMom));
@@ -245,7 +254,7 @@ class Integrator {
     MomentaField Mom2(P.Mom._grid);
     RealD RelativeError;
     Field diff(U._grid);
-    Real threshold = 1e-6;
+    Real threshold = 1e-8;
     int counter = 1;
     int MaxCounter = 100;
 
@@ -269,7 +278,7 @@ class Integrator {
       for (int mu = 0; mu < Nd; mu++) {
         auto Umu = PeekIndex<LorentzIndex>(U, mu);
         auto Pmu = PeekIndex<LorentzIndex>(sum, mu);
-        Umu = expMat(Pmu, ep * 0.5, 12) * Umu;
+        Umu = expMat(Pmu, ep * 0.5, 24) * Umu;
         PokeIndex<LorentzIndex>(NewU, ProjectOnGroup(Umu), mu);
       }
 
@@ -330,7 +339,7 @@ class Integrator {
   }
 
   void reverse_momenta(){
-    P.Mom *= 1.0;
+    P.Mom *= -1.0;
   }
 
   // to be used by the actionlevel class to iterate
diff --git a/tests/hmc/Test_hmc_WilsonGauge_Implicit.cc b/tests/hmc/Test_hmc_WilsonGauge_Implicit.cc
index 37a23f76..d9e04a32 100644
--- a/tests/hmc/Test_hmc_WilsonGauge_Implicit.cc
+++ b/tests/hmc/Test_hmc_WilsonGauge_Implicit.cc
@@ -53,14 +53,14 @@ int main(int argc, char **argv) {
   CheckpointerParameters CPparams;  
   CPparams.config_prefix = "ckpoint_lat";
   CPparams.rng_prefix = "ckpoint_rng";
-  CPparams.saveInterval = 20;
+  CPparams.saveInterval = 10;
   CPparams.format = "IEEE64BIG";
   
   TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
 
   RNGModuleParameters RNGpar;
   RNGpar.serial_seeds = "1 2 3 4 5";
-  RNGpar.parallel_seeds = "6 7 8 9 10 12";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
   TheHMC.Resources.SetRNGSeeds(RNGpar);
 
   // Construct observables
@@ -76,7 +76,7 @@ int main(int argc, char **argv) {
   // need wrappers of the fermionic classes 
   // that have a complex construction
   // standard
-  RealD beta = 0.0;
+  RealD beta = 5.6;
   WilsonGaugeActionR Waction(beta);
   
   ActionLevel<HMCWrapper::Field> Level1(1);
@@ -86,7 +86,7 @@ int main(int argc, char **argv) {
   /////////////////////////////////////////////////////////////
 
   // HMC parameters are serialisable 
-  TheHMC.Parameters.MD.MDsteps = 60;
+  TheHMC.Parameters.MD.MDsteps = 40;
   TheHMC.Parameters.MD.trajL   = 1.0;
 
   TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file