From 928fafd1f47a6449af466116951cbb011a56cc9c Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Mon, 15 Jan 2024 17:03:56 +0200
Subject: [PATCH] examples: compute diode effect properly in cpr_sns
Also try to clarify the mapping between alpha, eta and Gamma_DP/ST
---
examples/cpr_sns.py | 150 +++++++++++++++++++++++++++++++-------------
1 file changed, 105 insertions(+), 45 deletions(-)
diff --git a/examples/cpr_sns.py b/examples/cpr_sns.py
index 6bc1a6a..413356c 100644
--- a/examples/cpr_sns.py
+++ b/examples/cpr_sns.py
@@ -20,11 +20,11 @@ import collections
mem = joblib.Memory("cache")
-def get_solver(soc_alpha, eta, phi, L=10, h=0.5, D=1.0, n=5):
+def get_solver(soc_alpha, eta, phi, L=10, W=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.Ly = W
sol.Omega[...] = 0
sol.Delta[0, :] = np.eye(2) * np.exp(1j * phi / 2)
@@ -40,8 +40,8 @@ def get_solver(soc_alpha, eta, phi, L=10, h=0.5, D=1.0, n=5):
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)
+ Ax = soc_alpha * np.kron(S_0, S_y)
+ Ay = -soc_alpha * np.kron(S_0, S_x)
sol.Ux[...] = expm(1j * dx * Ax)
sol.Uy[...] = expm(1j * dy * Ay)
@@ -56,11 +56,10 @@ Res = collections.namedtuple("Res", ["x", "y", "J", "Jx", "Jy"])
@vectorize_parallel(returns_object=True, noarray=True)
@usadelndsoc.with_log_level(logging.WARNING)
-def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10):
+def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10, W=2):
+ sol = get_solver(soc_alpha=alpha_soc, eta=eta, L=L, W=W, n=n, phi=phi, h=h)
- sol = get_solver(soc_alpha=alpha_soc, eta=eta, L=L, n=n, phi=phi, h=h)
-
- E_typical = 10 + 10 * T
+ E_typical = 10 + 10 * abs(T)
w, a = get_matsubara_sum(T, E_typical)
t3 = np.kron(S_z, S_0)
@@ -76,51 +75,112 @@ def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10):
return (sol.x, sol.y, J, Jx, Jy)
-def main():
- T = 0.1
- alpha_soc = 0.9
- h = np.r_[0.0, 1.0]
- phi = np.linspace(-pi, pi, 37)
+def Gamma_to_alpha(Gamma_DP, Gamma_ST):
+ r"""
+ Transform from (Gamma_DP, Gamma_ST) to (alpha_soc, eta).
+
+ .. math::
+
+ \Gamma_{DP} &= 2 \alpha^2 p_F^2 \tau \\
+ \Gamma_{ST} &= \Gamma_{DP} \frac{\alpha\tau}{\xi_0}
- res = j(T, h[:, None], phi[None, :], alpha_soc=alpha_soc, L=3, n=20, mem=mem)
- res0 = j(T, h[:, None], phi[None, :], alpha_soc=alpha_soc, L=3, mem=mem)
- res1 = j(T, h[:, None], phi[None, :], alpha_soc=alpha_soc, L=3, n=10, mem=mem)
+ and :math:`\xi_0^2 = D/(2\pi\Delta_0)`.
+
+ The length unit is :math:`L_0 = \sqrt{D/E_0}` where
+ :math:`E_0` is the energy unit.
+ """
+
+ alpha_soc = (Gamma_DP / 4) ** 0.5
+ eta = Gamma_ST / (4 * np.sqrt(2 * pi) * alpha_soc**3)
+ return alpha_soc, eta
+
+
+def do(W_xi=6, multin=False):
+ T = 0.1
+ Gamma_DP = 10
+ Gamma_ST = 1
+ h = -10.0
+ phi = np.linspace(0, 2 * pi, 37)
+ xi = 1 / np.sqrt(2 * pi)
+ L = np.array([0.1, 1.5, 2, 3]) * xi
+ W = W_xi * xi
+
+ alpha, eta = Gamma_to_alpha(Gamma_DP, Gamma_ST)
+
+ res = j(
+ T, h, phi[None, :], alpha_soc=alpha, eta=eta, L=L[:, None], W=W, n=10, mem=mem
+ )
+
+ if multin:
+ res0 = j(
+ T,
+ h,
+ phi[None, :],
+ alpha_soc=alpha,
+ eta=eta,
+ L=L[:, None],
+ W=W,
+ n=5,
+ mem=mem,
+ )
+ res1 = j(
+ T,
+ h,
+ phi[None, :],
+ alpha_soc=alpha,
+ eta=eta,
+ L=L[:, None],
+ W=W,
+ n=20,
+ mem=mem,
+ )
Jx_mean = np.asarray(
[Res(*x).Jx[1:-1].sum(axis=1).mean(axis=0) for x in res.flat]
).reshape(res.shape)
- Jx0_mean = np.asarray(
- [Res(*x).Jx[1:-1].sum(axis=1).mean(axis=0) for x in res0.flat]
- ).reshape(res.shape)
+ if multin:
+ Jx0_mean = np.asarray(
+ [Res(*x).Jx[1:-1].sum(axis=1).mean(axis=0) for x in res0.flat]
+ ).reshape(res.shape)
- Jx1_mean = np.asarray(
- [Res(*x).Jx[1:-1].sum(axis=1).mean(axis=0) for x in res1.flat]
- ).reshape(res.shape)
+ Jx1_mean = np.asarray(
+ [Res(*x).Jx[1:-1].sum(axis=1).mean(axis=0) for x in res1.flat]
+ ).reshape(res.shape)
+
+ fig, axs = plt.subplots(1, 2, layout="compressed")
+
+ ax = axs[0]
+ ax.plot(phi / pi, Jx_mean.T / abs(Jx_mean).max(axis=1))
+ if multin:
+ ax.plot(phi / pi, Jx0_mean.T / abs(Jx_mean).max(axis=1), "k:")
+ ax.plot(phi / pi, Jx1_mean.T / abs(Jx_mean).max(axis=1), "k:")
+ ax.set_xlabel(r"$\varphi / \pi$")
+ ax.set_ylabel(r"$I / I_{\mathrm{max}}$")
+ ax.legend(L / xi, title=r"$L/\xi$", loc="lower right")
+
+ def eff(Jx):
+ Jm = Jx.min(axis=1)
+ Jp = Jx.max(axis=1)
+ return (abs(Jp) - abs(Jm)) / (abs(Jp) + abs(Jm))
- plt.plot(phi / pi, Jx_mean.T)
- plt.plot(phi / pi, Jx0_mean.T, 'k:')
- plt.plot(phi / pi, Jx1_mean.T, 'k:')
- plt.xlabel(r"$\varphi / \pi$")
- plt.ylabel(r"$I$")
- plt.legend(h, title="$h$")
- plt.savefig("cpr_sns.pdf")
-
- Jx_sym = (Jx_mean + Jx_mean[:,::-1])/2
- Jx0_sym = (Jx0_mean + Jx0_mean[:,::-1])/2
- Jx1_sym = (Jx1_mean + Jx1_mean[:,::-1])/2
- plt.clf()
- plt.plot(phi / pi, Jx_sym.T)
- plt.plot(phi / pi, Jx0_sym.T, 'k:')
- plt.plot(phi / pi, Jx1_sym.T, 'k:')
- plt.xlabel(r"$\varphi / \pi$")
- plt.ylabel(r"$[I(\varphi)+I(-\varphi)]/2$")
- plt.legend(h, title="$h$")
- plt.savefig("cpr_sns_sym.pdf")
-
- print(Jx_mean.mean(axis=1))
- print(Jx0_mean.mean(axis=1))
- print(Jx1_mean.mean(axis=1))
+ ax = axs[1]
+ ax.plot(L / xi, 100 * eff(Jx_mean))
+ if multin:
+ ax.plot(L / xi, 100 * eff(Jx0_mean), "k:")
+ ax.plot(L / xi, 100 * eff(Jx1_mean), "k:")
+ ax.set_xlabel(r"$L / \xi$")
+ ax.set_ylabel(r"$\eta$ [%]")
+
+ fig.suptitle(
+ rf"$W = {W/xi} \xi_0$ $\Gamma_{{DP}} = {Gamma_DP} \Delta_0$, $\Gamma_{{ST}} = {Gamma_ST} \Delta_0$ ($\tilde{{\eta}} = {eta:.3g}$, $\tilde{{\alpha}} = {alpha:.3g}$)"
+ )
+
+ fig.savefig("cpr_sns.pdf")
+
+
+def main():
+ do()
if __name__ == "__main__":
--
GitLab