diff --git a/tests/test_basic.py b/tests/test_basic.py
index cf7cffae8523a433fc01e45835ada1445127837b..caded57da05037739cfa9b93ef07e5e4cce4c49b 100644
--- a/tests/test_basic.py
+++ b/tests/test_basic.py
@@ -201,3 +201,12 @@ def test_gauge_invariance(omega, use_numdifftools):
         D = numdifftools.Derivative(S_s)
         dS_ds = D(0.0) / 16
         assert_allclose(div_J, dS_ds, rtol=1e-12)
+
+    # The derivative must also match the gradient
+    dS_dPhi = s._core.grad(Phi).reshape(Phi.shape)
+    dPhi = 0 * Phi
+    dPhi[i, j, 0] = 2 * Phi[i, j, 0]
+    dPhi[i, j, 1] = -2 * Phi[i, j, 1]
+    dS_ds2 = (dS_dPhi * dPhi).sum() / 16
+
+    assert_allclose(div_J, dS_ds2, rtol=1e-12)