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 d25f324f..a6d401a1 100644
--- a/src/mp00ac.f90
+++ b/src/mp00ac.f90
@@ -532,28 +532,34 @@ subroutine mp00ac( Ndof, Xdof, Fdof, Ddof, Ldfjac, iflag ) ! argument list is fi
 
 
   ! if (GaugeInvariantH) then
+  dH1 = zero
+  dH2 = zero
   do ll = 0, Lrad(lvol)
     !             ideriv=0   mn=0
     id = Ate(lvol,0        , 1)%i(ll)
-    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)
+    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 (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)
+    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