diff --git a/lib/qcd/action/gauge/GaugeImpl.h b/lib/qcd/action/gauge/GaugeImpl.h
index 8e9ad36f..fd3bf55a 100644
--- a/lib/qcd/action/gauge/GaugeImpl.h
+++ b/lib/qcd/action/gauge/GaugeImpl.h
@@ -96,6 +96,17 @@ public:
 
   }
 
+  static inline RealD FieldSquareNorm(Field& U){
+    LatticeComplex Hloc(U._grid);
+    Hloc = zero;
+    for (int mu = 0; mu < Nd; mu++) {
+      auto Umu = PeekIndex<LorentzIndex>(U, mu);
+      Hloc += trace(Umu * Umu);
+    }
+    Complex Hsum = sum(Hloc);
+    return Hsum.real();
+  }
+
 };
 
 // Composition with smeared link, bc's etc.. probably need multiple inheritance
diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index 09e1c275..d2c1aef1 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -253,17 +253,7 @@ class Integrator {
   // Calculate action
   RealD S(Field& U) {  // here also U not used
 
-    LatticeComplex Hloc(U._grid);
-    Hloc = zero;
-    // Momenta  --- modify this (take the trace of field)
-    // momenta_action()
-    for (int mu = 0; mu < Nd; mu++) {
-      auto Pmu = PeekIndex<LorentzIndex>(P, mu);
-      Hloc -= trace(Pmu * Pmu);
-    }
-    Complex Hsum = sum(Hloc);
-
-    RealD H = Hsum.real();
+    RealD H = - FieldImplementation::FieldSquareNorm(P); // - trace (P*P)
     RealD Hterm;
     std::cout << GridLogMessage << "Momentum action H_p = " << H << "\n";