Skip to content
Snippets Groups Projects
Commit da76e299 authored by patavirt's avatar patavirt
Browse files

tests: fixup prefactors

parent f5412af5
No related branches found
No related tags found
No related merge requests found
......@@ -270,7 +270,7 @@ def test_soc_analytic(n=10):
sz = np.array([[1, 0], [0, -1]])
s0 = np.array([[1, 0], [0, 1]])
def get_sol(soc_alpha, alpha, L=10, h=0.5):
def get_sol(soc_alpha, alpha, L=10, h=0.5, D=1.0):
sol = Solver(nx=n, ny=n)
sol.mask[...] = MASK_NONE
sol.Lx = L
......@@ -279,6 +279,7 @@ def test_soc_analytic(n=10):
sol.Omega[...] = 0
sol.Delta[:, :] = np.eye(2)
sol.D = D
sol.alpha = alpha
nx, ny = sol.shape
......@@ -295,10 +296,16 @@ def test_soc_analytic(n=10):
return sol
sol = get_sol(soc_alpha=1.0, alpha=1.0)
sol0 = get_sol(soc_alpha=0, alpha=0)
# Solver without SOC
h = 0.5
sol0 = get_sol(soc_alpha=0, alpha=0, h=h)
T = 0.2
# 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)
......@@ -308,18 +315,35 @@ def test_soc_analytic(n=10):
Fs = []
Gs = []
Js = []
Jans = []
dd = (sol.Lx / sol.shape[0], sol.Ly / sol.shape[1])
Jtot = 0
Jantot = 0
for ww, aa in zip(w, a):
res0 = sol0.solve(omega=ww, finite_difference=False)
res0 = sol0.solve(omega=ww)
sol.omega = ww
J = sol._core.grad_A(res0.Phi.copy()).transpose(0, 1, 2, 4, 3) * (-1 / 8)
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
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)
)
Jantot += aa * Jy_an
Fs.append(res0.F.copy())
Gs.append(res0.G.copy())
Phis.append(res0.Phi.copy())
Js.append(J)
Jans.append(Jy_an)
return locals()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment