From 42dab0ae7a875593f97d17a7c2d52a1b1c5a5c75 Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Wed, 21 Sep 2022 14:59:36 +0300
Subject: [PATCH] examples: add initial cpr_sns
---
examples/cpr_sns.py | 96 +++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 96 insertions(+)
create mode 100644 examples/cpr_sns.py
diff --git a/examples/cpr_sns.py b/examples/cpr_sns.py
new file mode 100644
index 0000000..371a2ce
--- /dev/null
+++ b/examples/cpr_sns.py
@@ -0,0 +1,96 @@
+# -*- 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()
--
GitLab