tests: make test_soc_analytic a test

tests/test_solver.py

@@ -264,7 +264,10 @@ def curve_separation(x1, y1, x2, y2):
     return s.min(axis=1).max()
 
 
-def test_soc_analytic(n=10):
+@pytest.mark.parametrize(
+    "T,n,alpha_soc", [(0.1, 10, 0.01), (0.4, 10, 0.01), (0.4, 5, 0.0123)]
+)
+def test_soc_analytic(T, n, alpha_soc):
     sx = np.array([[0, 1], [1, 0]])
     sy = np.array([[0, -1j], [1j, 0]])
     sz = np.array([[1, 0], [0, -1]])
@@ -301,11 +304,9 @@ def test_soc_analytic(n=10):
     sol0 = get_sol(soc_alpha=0, alpha=0, h=h)
 
     # Solver for evaluating current perturbatively in SOC
-    alpha_soc = 0.01
     alpha = 1.0
     sol = get_sol(soc_alpha=alpha_soc, alpha=alpha, h=h, D=0)
 
-    T = 0.4
     E_typical = 10 + 10 * T
     w, a = get_matsubara_sum(T, E_typical)
 
@@ -325,20 +326,24 @@ def test_soc_analytic(n=10):
         res0 = sol0.solve(omega=ww)
         sol.omega = ww
         J = sol._core.grad_A(res0.Phi.copy()).transpose(0, 1, 2, 4, 3) * (1 / 16)
-        J[:, :, 0] /= dd[1] * sol.Ly
-        J[:, :, 1] /= dd[0] * sol.Lx
-        J[:, :, 2] /= dd[1] * sol.Ly
-        J[:, :, 3] /= dd[0] * sol.Lx
-        Jtot += aa * J
+        J[:, :, 0] /= dd[1]
+        J[:, :, 1] /= dd[0]
+        J[:, :, 2] /= dd[1]
+        J[:, :, 3] /= dd[0]
+        Jtot += pi * aa * J
 
         gp = (ww - 1j * h) / np.sqrt((ww - 1j * h) ** 2 + 1)
         gm = (ww + 1j * h) / np.sqrt((ww + 1j * h) ** 2 + 1)
         fp = 1 / np.sqrt((ww - 1j * h) ** 2 + 1)
         fm = 1 / np.sqrt((ww + 1j * h) ** 2 + 1)
         Jy_an = (
-            (gp - gm) * (1 + gp * gm + fp * fm) * (-1j * pi * alpha * alpha_soc**3)
+            (gp - gm) * (1 + gp * gm + fp * fm) * (-1j * alpha_soc**3) * (32 * alpha)
         )
-        Jantot += aa * Jy_an
+        Jantot += pi * aa * Jy_an
+
+        Jy = tr(J @ t3)[n // 2, n // 2, 1]
+
+        assert_allclose(-Jy, Jy_an, rtol=1e-2)
 
         Fs.append(res0.F.copy())
         Gs.append(res0.G.copy())