add files

usadelndsoc/bcs.py

+# -*- mode:python; coding: utf-8; eval: (blacken-mode) -*-
+"""
+BCS superconductor gap
+"""
+# Author: Pauli Virtanen
+# License: GNU Affero General Public License, see LICENSE.txt for details
+import numpy as np
+from numpy import pi, exp, sqrt, log, euler_gamma
+from scipy.optimize import brenth, golden
+from scipy.interpolate import make_interp_spline, PPoly
+
+from . import matsubara
+
+__all__ = ["BCS_Delta"]
+
+
+@np.vectorize
+def BCS_Delta(T, Tc, h=0):
+    """
+    Evaluate Delta in a bulk BCS superconductor.
+    """
+    Delta_0 = pi / exp(euler_gamma) * Tc
+
+    if T < 1e-2 * Delta_0:
+        T = 1e-2 * Delta_0
+
+    w, a = matsubara.get_matsubara_sum(T, 100 * max(Tc, abs(h)), steps=128, max_ne=3000)
+
+    def delta_fun(Delta):
+        Omega = w + 1j * h
+        r = 1 / abs(w) - (1 / sqrt(Omega**2 + abs(Delta) ** 2)).real
+        res = log(Tc / T) - pi * (a * r).sum(axis=-1)
+        return res
+
+    def F_fun(Delta):
+        Omega = w + 1j * h
+        r = np.real(
+            2 * Omega * np.sign(Omega.real)
+            - (2 * Omega**2 + Delta**2) / np.sqrt(Omega**2 + Delta**2)
+        )
+        return pi * (a * r).sum(axis=-1)
+
+    x = np.linspace(0, 2 * Delta_0, 500)
+
+    # Bracket roots
+    y = delta_fun(x[:, None])
+    spl = make_interp_spline(x, y, k=1)
+    spl = PPoly.from_spline(spl, extrapolate=False)
+    r = spl.roots()
+
+    if r.size == 0:
+        return 0.0
+
+    r_a = np.r_[0, (r[:-1] + r[1:]) / 2]
+    r_b = np.r_[r_a[1:], r[-1] + 0.1]
+
+    # Refine roots and find the one with smallest free energy
+    r_min = 0
+    F_min = F_fun(0.0)
+
+    for aa, bb in zip(r_a, r_b):
+        try:
+            r = brenth(delta_fun, aa, bb)
+        except ValueError:
+            print("error:", T, Tc, h, aa, bb, delta_fun(aa), delta_fun(bb))
+            raise
+        F = F_fun(r)
+        if F < F_min:
+            r_min = r
+
+    return r_min

usadelndsoc/meson.build

@@ -7,7 +7,7 @@ py3.extension_module(
 )
 
 py3.install_sources(
-  ['__init__.py', 'matsubara.py'],
+  ['__init__.py', 'bcs.py', 'matsubara.py', 'solver.py'],
   pure: false,
   subdir : 'usadelndsoc',
 )