From dc36d272ce1701cd4c9e8bb00b1509bc03008c37 Mon Sep 17 00:00:00 2001
From: Chulwoo Jung <chulwoo@bnl.gov>
Date: Tue, 21 Nov 2023 13:48:51 -0500
Subject: [PATCH] Gauge RMHMC conserving dH

---
 Grid/qcd/action/gauge/GaugeImplTypes.h    |  2 +-
 Grid/qcd/action/gauge/WilsonGaugeAction.h |  2 +-
 Grid/qcd/action/scalar/ScalarImpl.h       |  2 +-
 Grid/qcd/hmc/HMC.h                        |  2 +-
 Grid/qcd/hmc/integrators/Integrator.h     | 40 +++++++++++++++++++----
 Grid/qcd/utils/Metric.h                   |  3 +-
 HMC/Mobius2p1p1fEOFA_4Gev.cc              |  2 +-
 7 files changed, 41 insertions(+), 12 deletions(-)

diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h
index a9af1fae..3aafc57e 100644
--- a/Grid/qcd/action/gauge/GaugeImplTypes.h
+++ b/Grid/qcd/action/gauge/GaugeImplTypes.h
@@ -32,7 +32,7 @@ directory
 
 NAMESPACE_BEGIN(Grid);
 
-#define CPS_MD_TIME
+#undef CPS_MD_TIME
 
 #ifdef CPS_MD_TIME
 #define HMC_MOMENTUM_DENOMINATOR (2.0)
diff --git a/Grid/qcd/action/gauge/WilsonGaugeAction.h b/Grid/qcd/action/gauge/WilsonGaugeAction.h
index 0156811d..ca1e27bc 100644
--- a/Grid/qcd/action/gauge/WilsonGaugeAction.h
+++ b/Grid/qcd/action/gauge/WilsonGaugeAction.h
@@ -73,7 +73,7 @@ public:
 	////////// boundary phase /////////////
       auto pha = Params.boundary_phases[mu];
       scalar_type phase( real(pha),imag(pha) );
-   std::cout<< GridLogMessage << "[WilsonGaugeAction] boundary "<<mu<<" "<<phase<< std::endl; 
+      std::cout<< GridLogIterative << "[WilsonGaugeAction] boundary "<<mu<<" "<<phase<< std::endl; 
 
 	int L   = GaugeGrid->GlobalDimensions()[mu];
         int Lmu = L - 1;
diff --git a/Grid/qcd/action/scalar/ScalarImpl.h b/Grid/qcd/action/scalar/ScalarImpl.h
index 13bd6c90..dafc520b 100644
--- a/Grid/qcd/action/scalar/ScalarImpl.h
+++ b/Grid/qcd/action/scalar/ScalarImpl.h
@@ -1,6 +1,6 @@
 #pragma once
 
-#define CPS_MD_TIME 
+#undef CPS_MD_TIME 
 
 #ifdef CPS_MD_TIME
 #define HMC_MOMENTUM_DENOMINATOR (2.0)
diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h
index f91f03e3..4151e1c2 100644
--- a/Grid/qcd/hmc/HMC.h
+++ b/Grid/qcd/hmc/HMC.h
@@ -209,7 +209,7 @@ private:
     // updated state action
     //////////////////////////////////////////////////////////////////////////////////////////////////////
     std::cout << GridLogMessage << "--------------------------------------------------\n";
-    std::cout << GridLogMessage << "Compute final action";
+    std::cout << GridLogMessage << "Compute final action" <<std::endl;
     RealD H1 = TheIntegrator.S(U);  
     std::cout << GridLogMessage << "--------------------------------------------------\n";
 
diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h
index 1673905d..292f001e 100644
--- a/Grid/qcd/hmc/integrators/Integrator.h
+++ b/Grid/qcd/hmc/integrators/Integrator.h
@@ -91,6 +91,8 @@ public:
   RepresentationPolicy Representations;
   IntegratorParameters Params;
 
+  RealD Saux,Smom,Sg;
+
   //Filters allow the user to manipulate the conjugate momentum, for example to freeze links in DDHMC
   //It is applied whenever the momentum is updated / refreshed
   //The default filter does nothing
@@ -384,7 +386,8 @@ public:
       P(grid, M),
       levels(Aset.size()),
       Smearer(Sm),
-      Representations(grid) 
+      Representations(grid),
+      Saux(0.),Smom(0.),Sg(0.)
   {
     t_P.resize(levels, 0.0);
     t_U = 0.0;
@@ -524,7 +527,10 @@ public:
     std::cout << GridLogIntegrator << "Integrator refresh" << std::endl;
 
     std::cout << GridLogIntegrator << "Generating momentum" << std::endl;
-    FieldImplementation::generate_momenta(P.Mom, sRNG, pRNG);
+//    FieldImplementation::generate_momenta(P.Mom, sRNG, pRNG);
+    P.M.ImportGauge(U);
+    P.MomentaDistribution(sRNG,pRNG);
+
 
     // Update the smeared fields, can be implemented as observer
     // necessary to keep the fields updated even after a reject
@@ -579,9 +585,22 @@ public:
 
     std::cout << GridLogIntegrator << "Integrator action\n";
 
-    RealD H = - FieldImplementation::FieldSquareNorm(P.Mom)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+//    RealD H = - FieldImplementation::FieldSquareNorm(P.Mom)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+//    RealD Hterm;
+
+//    static RealD Saux=0.,Smom=0.,Sg=0.;
+
+    RealD H = - FieldImplementation::FieldSquareNorm(P.Mom)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+    std::cout << GridLogMessage << "S:FieldSquareNorm H_p = " << H << "\n";
+    std::cout << GridLogMessage << "S:dSField = " << H-Smom << "\n";
+    Smom=H;
+    P.M.ImportGauge(U);
+    RealD Hterm = - P.MomentaAction();
+    std::cout << GridLogMessage << "S:Momentum action H_p = " << Hterm << "\n";
+    std::cout << GridLogMessage << "S:dSMom = " << Hterm-Saux << "\n";
+    Saux=Hterm;
+    H = Hterm;
 
-    RealD Hterm;
 
     // Actions
     for (int level = 0; level < as.size(); ++level) {
@@ -623,9 +642,18 @@ public:
 
     std::cout << GridLogIntegrator << "Integrator initial action\n";
 
+//    RealD H = - FieldImplementation::FieldSquareNorm(P.Mom)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+//    RealD Hterm;
     RealD H = - FieldImplementation::FieldSquareNorm(P.Mom)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
-
-    RealD Hterm;
+    std::cout << GridLogMessage << "S:FieldSquareNorm H_p = " << H << "\n";
+    std::cout << GridLogMessage << "S:dSField = " << H-Smom << "\n";
+    Smom=H;
+    P.M.ImportGauge(U);
+    RealD Hterm = - P.MomentaAction();
+    std::cout << GridLogMessage << "S:Momentum action H_p = " << Hterm << "\n";
+    std::cout << GridLogMessage << "S:dSMom = " << Hterm-Saux << "\n";
+    Saux=Hterm;
+    H = Hterm;
 
     // Actions
     for (int level = 0; level < as.size(); ++level) {
diff --git a/Grid/qcd/utils/Metric.h b/Grid/qcd/utils/Metric.h
index f09cb3fc..c063aadd 100644
--- a/Grid/qcd/utils/Metric.h
+++ b/Grid/qcd/utils/Metric.h
@@ -232,7 +232,8 @@ public:
   }
 
   void update_auxiliary_momenta(RealD ep){
-//    if(!M.Trivial()) 
+      std::cout << GridLogIntegrator << "AuxMom update_auxiliary_fields: " << std::sqrt(norm2(AuxMom)) << std::endl;
+      std::cout << GridLogIntegrator << "AuxField update_auxiliary_fields: " << std::sqrt(norm2(AuxField)) << std::endl;
     {
       AuxMom -= ep * AuxField * HMC_MOMENTUM_DENOMINATOR;
       std::cout << GridLogIntegrator << "AuxMom update_auxiliary_fields: " << std::sqrt(norm2(AuxMom)) << std::endl;
diff --git a/HMC/Mobius2p1p1fEOFA_4Gev.cc b/HMC/Mobius2p1p1fEOFA_4Gev.cc
index 02065a0c..b5aa8d8d 100644
--- a/HMC/Mobius2p1p1fEOFA_4Gev.cc
+++ b/HMC/Mobius2p1p1fEOFA_4Gev.cc
@@ -299,7 +299,7 @@ int main(int argc, char **argv) {
   
   double ActionStoppingCondition     = 1e-8;
   double DerivativeStoppingCondition = 1e-6;
-  double MaxCGIterations = 30000;
+  double MaxCGIterations =  100000;
 
   ////////////////////////////////////
   // Collect actions