diff --git a/ci/G1V03L2Fi/compare.h5 b/ci/G1V03L2Fi/compare.h5
index 10297c33..f7579e2a 100644
Binary files a/ci/G1V03L2Fi/compare.h5 and b/ci/G1V03L2Fi/compare.h5 differ
diff --git a/ci/G2V32L1Fi/compare.h5 b/ci/G2V32L1Fi/compare.h5
index 3c41ae2a..2169e49d 100644
Binary files a/ci/G2V32L1Fi/compare.h5 and b/ci/G2V32L1Fi/compare.h5 differ
diff --git a/ci/G3V08L3Fi/compare.h5 b/ci/G3V08L3Fi/compare.h5
index 53c02cd0..0afc2c3d 100644
Binary files a/ci/G3V08L3Fi/compare.h5 and b/ci/G3V08L3Fi/compare.h5 differ
diff --git a/ci/G3V08L3Fr/compare.h5 b/ci/G3V08L3Fr/compare.h5
index b6664412..21782275 100644
Binary files a/ci/G3V08L3Fr/compare.h5 and b/ci/G3V08L3Fr/compare.h5 differ
diff --git a/ci/toroidal_freeboundary_vacuum/compare.h5 b/ci/toroidal_freeboundary_vacuum/compare.h5
index 430a6e0f..39c1692f 100644
Binary files a/ci/toroidal_freeboundary_vacuum/compare.h5 and b/ci/toroidal_freeboundary_vacuum/compare.h5 differ
diff --git a/src/mp00ac.f90 b/src/mp00ac.f90
index 2ecaf1ab..a6d401a1 100644
--- a/src/mp00ac.f90
+++ b/src/mp00ac.f90
@@ -117,7 +117,7 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
-  use constants, only : zero, half, one
+  use constants, only : zero, half, one, pi2
 
   use numerical, only : small, machprec
 
@@ -127,7 +127,8 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi
                         mu, helicity, iota, oita, curtor, curpol, Lrad, Ntor,&
                        !Lposdef, &
                         Lconstraint, mupftol, &
-                        Lmatsolver, NiterGMRES, epsGMRES, LGMRESprec, epsILU
+                        Lmatsolver, NiterGMRES, epsGMRES, LGMRESprec, epsILU,&
+                        mpol
 
   use cputiming, only : Tmp00ac
 
@@ -196,6 +197,11 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi
   REAL, allocatable    :: wk(:)
   INTEGER,allocatable  :: jw(:), iperm(:)
 
+  REAL               :: cheby(0:Lrad(ivol),0:2), zernike(0:Lrad(1),0:mpol,0:2)
+  REAL               :: dH1, dH2
+  INTEGER            :: ll, id
+
+
   BEGIN(mp00ac)
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
@@ -524,6 +530,39 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi
     lABintegral(lvol) = half * sum( solution(1:NN,0) * Ddotx(1:NN) )
   endif
 
+
+  ! if (GaugeInvariantH) then
+  dH1 = zero
+  dH2 = zero
+  do ll = 0, Lrad(lvol)
+    !             ideriv=0   mn=0
+    id = Ate(lvol,0        , 1)%i(ll)
+    if (id/=0) then ! Indirection maps some elements to 0, to "discard" them during matrix construction. Ignore them here.
+      if (Lcoordinatesingularity) then
+        call get_zernike_d2(small, Lrad(lvol), mpol, zernike)
+        dH1 = dH1 + pi2 * solution(id,0) * zernike(ll,0,0)
+      else
+        call get_cheby_d2(-one+small, Lrad(lvol), cheby(0:Lrad(lvol),0:2))
+        dH1 = dH1 + pi2 * solution(id,0) * cheby(ll,0)
+      endif
+    endif
+    
+    !             ideriv=0   mn=0
+    id = Aze(lvol,0        , 1)%i(ll)
+    if (id/=0) then ! Indirection maps some elements to 0, to "discard" them during matrix construction. Ignore them here.
+      if (Lcoordinatesingularity) then
+        call get_zernike_d2(one, Lrad(lvol), mpol, zernike)
+        dH2 = dH2 + pi2 * solution(id,0) * zernike(ll,0,0)
+      else
+        call get_cheby_d2(one, Lrad(lvol), cheby(0:Lrad(lvol),0:2))
+        dH2 = dH2 + pi2 * solution(id,0) * cheby(ll,0)
+      endif
+    endif
+  enddo
+  
+  lABintegral(lvol) = lABintegral(lvol) - dH1 - dH2
+  ! endif
+
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
   DALLOCATE( matrix )
   DALLOCATE( rhs    )