examples: add initial cpr_sns

examples/cpr_sns.py

+# -*- coding:utf-8; eval: (blacken-mode) -*-
+"""
+Calculate anomalous current in bulk case
+"""
+import joblib
+import numpy as np
+from numpy import pi
+import matplotlib.pyplot as plt
+from scipy.linalg import expm
+import logging
+import contextlib
+
+import usadelndsoc
+from usadelndsoc.matsubara import get_matsubara_sum
+from usadelndsoc.solver import *
+from usadelndsoc.util import vectorize_parallel
+
+mem = joblib.Memory("cache")
+
+
+def get_solver(soc_alpha, eta, phi, L=10, h=0.5, D=1.0, n=5):
+    sol = Solver(nx=n, ny=n)
+    sol.mask[...] = MASK_NONE
+    sol.Lx = L
+    sol.Ly = L
+
+    sol.Omega[...] = 0
+    sol.Delta[0, :] = np.eye(2) * np.exp(1j * phi / 2)
+    sol.Delta[-1, :] = np.eye(2) * np.exp(-1j * phi / 2)
+
+    sol.mask[0, :] = MASK_TERMINAL
+    sol.mask[-1, :] = MASK_TERMINAL
+
+    sol.D = D
+    sol.eta = eta
+
+    nx, ny = sol.shape
+    dx = sol.Lx / nx
+    dy = sol.Ly / ny
+
+    Ax = 2 * soc_alpha * np.kron(S_0, S_y)
+    Ay = -2 * soc_alpha * np.kron(S_0, S_x)
+
+    sol.Ux[...] = expm(1j * dx * Ax)
+    sol.Uy[...] = expm(1j * dy * Ay)
+
+    sol.Omega[...] += h * np.kron(S_z, S_y)
+
+    return sol
+
+
+# @usadelndsoc.with_log_level(logging.WARNING)
+@vectorize_parallel(returns_object=True, noarray=True)
+def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10, perturbative=False):
+    sol = get_solver(soc_alpha=alpha_soc, eta=eta, L=L, n=n, phi=phi, h=h)
+
+    E_typical = 10 + 10 * T
+    w, a = get_matsubara_sum(T, E_typical)
+
+    t3 = np.kron(S_z, S_0)
+    res = sol.solve_many(omega=w)
+
+    t3 = np.diag([1, 1, -1, -1])
+    J = tr(res.J @ t3) * (-1 / 16)
+    J = -1j * pi * (J * a[:, None, None, None]).sum(axis=0)
+
+    Jx = (J[:, :, 0] + J[:, :, 2]).real / 2
+    Jy = (J[:, :, 1] + J[:, :, 3]).real / 2
+
+    return sol.x, sol.y, J, Jx, Jy
+
+
+def main():
+    T = 0.02
+    alpha_soc = 0.05
+    Da2 = alpha_soc**2
+    h = np.linspace(0, 1.75, 29)
+
+    res = j_anom(T, h, alpha_soc=alpha_soc, perturbative=False, mem=mem)
+    j = np.asarray([x[0] for x in res])
+    j_an = np.asarray([x[1] for x in res])
+
+    t3 = np.diag([1, 1, -1, -1])
+    j = tr(j @ t3)[:, 2, 2, UP]
+
+    plt.plot(h, -j.real, "-", label=r"numerics")
+    plt.plot(h, -j_an.real, "--", label=r"analytic")
+    plt.title(rf"$T = {T} \Delta$, $D\alpha^2 = {Da2:.4} \Delta$")
+    plt.xlabel(r"$h/\Delta$")
+    plt.ylabel(r"$j_{\mathrm{an}}$")
+    plt.legend()
+    plt.savefig("bulk_j_an.pdf")
+
+
+if __name__ == "__main__":
+    main()