Gauge RMHMC conserving dH

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)

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;

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)

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";
 

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
-
+    std::cout << GridLogMessage << "S:FieldSquareNorm H_p = " << H << "\n";
-    RealD Hterm;
+    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) {

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;

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