diff --git a/desc/integrals/singularities.py b/desc/integrals/singularities.py
index 52702ba2f..c389cc439 100644
--- a/desc/integrals/singularities.py
+++ b/desc/integrals/singularities.py
@@ -1117,14 +1117,7 @@ def LHS(Phi_mn):
     return Phi_mn, Phi_transform
 
 
-def compute_grad_Phi(
-    eq,
-    eval_grid,
-    source_grid,
-    Phi_mn,
-    basis,
-    return_data=False,
-):
+def compute_grad_Phi(eq, eval_grid, source_grid, Phi_mn, basis):
     """Computes vacuum field ∇Φ on ∂D.
 
     Let D, D^∁ denote the interior, exterior of a toroidal region with
@@ -1206,8 +1199,6 @@ def compute_grad_Phi(
             loop=True,
         )
     )
-    if return_data:
-        return grad_Phi, evl_data
     return grad_Phi
 
 
diff --git a/tests/test_integrals.py b/tests/test_integrals.py
index 6f7068040..4c4b27bf6 100644
--- a/tests/test_integrals.py
+++ b/tests/test_integrals.py
@@ -743,7 +743,7 @@ def test(G):
                 eq=eq, B0n=B0n, source_grid=src_grid, Phi_N=max(eq.N, 1)
             )
             evl_grid = Phi_transform.grid
-            grad_Phi, evl_data = compute_grad_Phi(
+            grad_Phi = compute_grad_Phi(
                 eq=eq,
                 eval_grid=evl_grid,
                 source_grid=src_grid,
@@ -751,7 +751,8 @@ def test(G):
                 basis=Phi_transform.basis,
                 return_data=True,
             )
-            dPhi_dn = dot(grad_Phi, evl_data["n_rho"])
+            n_rho = eq.compute("n_rho", grid=evl_grid)["n_rho"]
+            dPhi_dn = dot(grad_Phi, n_rho)
             B0n, _ = B0.compute_Bnormal(
                 eq.surface, eval_grid=evl_grid, source_grid=src_grid
             )