diff --git a/phase_diagram.ipynb b/phase_diagram.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..4b73362b194d2c8f5df7773ac40ef0c9ef55e929
--- /dev/null
+++ b/phase_diagram.ipynb
@@ -0,0 +1,616 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%%writefile phase_diagram.py\n",
+ "\n",
+ "\n",
+ "import inspect\n",
+ "import operator\n",
+ "import sys\n",
+ "from collections import OrderedDict\n",
+ "from copy import copy, deepcopy\n",
+ "from functools import wraps\n",
+ "from itertools import product\n",
+ "from types import SimpleNamespace\n",
+ "\n",
+ "import kwant\n",
+ "import numpy as np\n",
+ "import scipy.constants\n",
+ "from kwant.continuum.discretizer import discretize\n",
+ "from skimage import measure\n",
+ "\n",
+ "import adaptive\n",
+ "import pfaffian as pf\n",
+ "\n",
+ "assert sys.version_info >= (3, 6), \"Use Python ≥3.6\"\n",
+ "\n",
+ "# Parameters taken from arXiv:1204.2792\n",
+ "# All constant parameters, mostly fundamental constants, in a SimpleNamespace.\n",
+ "constants = SimpleNamespace(\n",
+ " m_eff=0.015 * scipy.constants.m_e, # effective mass in kg\n",
+ " hbar=scipy.constants.hbar,\n",
+ " m_e=scipy.constants.m_e,\n",
+ " eV=scipy.constants.eV,\n",
+ " e=scipy.constants.e,\n",
+ " c=1e18 / (scipy.constants.eV * 1e-3), # to get to meV * nm^2\n",
+ " mu_B=scipy.constants.physical_constants[\"Bohr magneton in eV/T\"][0] * 1e3,\n",
+ ")\n",
+ "\n",
+ "constants.t = (constants.hbar ** 2 / (2 * constants.m_eff)) * constants.c\n",
+ "\n",
+ "\n",
+ "def get_names(sig):\n",
+ " names = [\n",
+ " (name, value)\n",
+ " for name, value in sig.parameters.items()\n",
+ " if value.kind\n",
+ " in (inspect.Parameter.POSITIONAL_OR_KEYWORD, inspect.Parameter.KEYWORD_ONLY)\n",
+ " ]\n",
+ " return OrderedDict(names)\n",
+ "\n",
+ "\n",
+ "def filter_kwargs(sig, names, kwargs):\n",
+ " names_in_kwargs = [(name, value) for name, value in kwargs.items() if name in names]\n",
+ " return OrderedDict(names_in_kwargs)\n",
+ "\n",
+ "\n",
+ "def skip_pars(names1, names2, num_skipped):\n",
+ " skipped_pars1 = list(names1.keys())[:num_skipped]\n",
+ " skipped_pars2 = list(names2.keys())[:num_skipped]\n",
+ " if skipped_pars1 == skipped_pars2:\n",
+ " pars1 = list(names1.values())[num_skipped:]\n",
+ " pars2 = list(names2.values())[num_skipped:]\n",
+ " else:\n",
+ " raise Exception(\"First {} arguments \" \"have to be the same\".format(num_skipped))\n",
+ " return pars1, pars2\n",
+ "\n",
+ "\n",
+ "def combine(f, g, operator, num_skipped=0):\n",
+ " if not callable(f) or not callable(g):\n",
+ " raise Exception(\"One of the functions is not a function\")\n",
+ "\n",
+ " sig1 = inspect.signature(f)\n",
+ " sig2 = inspect.signature(g)\n",
+ "\n",
+ " names1 = get_names(sig1)\n",
+ " names2 = get_names(sig2)\n",
+ "\n",
+ " pars1, pars2 = skip_pars(names1, names2, num_skipped)\n",
+ " skipped_pars = list(names1.values())[:num_skipped]\n",
+ "\n",
+ " pars1_names = {p.name for p in pars1}\n",
+ " pars2 = [p for p in pars2 if p.name not in pars1_names]\n",
+ "\n",
+ " parameters = pars1 + pars2\n",
+ " kind = inspect.Parameter.POSITIONAL_OR_KEYWORD\n",
+ " parameters = [p.replace(kind=kind) for p in parameters]\n",
+ " parameters = skipped_pars + parameters\n",
+ "\n",
+ " def wrapped(*args):\n",
+ " d = {p.name: arg for arg, p in zip(args, parameters)}\n",
+ " fval = f(*[d[name] for name in names1.keys()])\n",
+ " gval = g(*[d[name] for name in names2.keys()])\n",
+ " return operator(fval, gval)\n",
+ "\n",
+ " wrapped.__signature__ = inspect.Signature(parameters=parameters)\n",
+ " return wrapped\n",
+ "\n",
+ "\n",
+ "def memoize(obj):\n",
+ " cache = obj.cache = {}\n",
+ "\n",
+ " @wraps(obj)\n",
+ " def memoizer(*args, **kwargs):\n",
+ " key = str(args) + str(kwargs)\n",
+ " if key not in cache:\n",
+ " cache[key] = obj(*args, **kwargs)\n",
+ " return cache[key]\n",
+ "\n",
+ " return memoizer\n",
+ "\n",
+ "\n",
+ "@memoize\n",
+ "def discretized_hamiltonian(a, which_lead=None):\n",
+ " ham = (\n",
+ " \"(0.5 * hbar**2 * (k_x**2 + k_y**2 + k_z**2) / m_eff * c - mu + V) * kron(sigma_0, sigma_z) + \"\n",
+ " \"alpha * (k_y * kron(sigma_x, sigma_z) - k_x * kron(sigma_y, sigma_z)) + \"\n",
+ " \"0.5 * g * mu_B * (B_x * kron(sigma_x, sigma_0) + B_y * kron(sigma_y, sigma_0) + B_z * kron(sigma_z, sigma_0)) + \"\n",
+ " \"Delta * kron(sigma_0, sigma_x)\"\n",
+ " )\n",
+ "\n",
+ " subst_sm = {\"Delta\": 0}\n",
+ "\n",
+ " if which_lead is not None:\n",
+ " subst_sm[\"V\"] = f\"V_{which_lead}(z, V_0, V_r, V_l, x0, sigma, r1)\"\n",
+ " subst_sm[\"mu\"] = f\"mu_{which_lead}(x0, sigma, mu_lead, mu_wire)\"\n",
+ " else:\n",
+ " subst_sm[\"V\"] = \"V(x, z, V_0, V_r, V_l, x0, sigma, r1)\"\n",
+ " subst_sm[\"mu\"] = \"mu(x, x0, sigma, mu_lead, mu_wire)\"\n",
+ "\n",
+ " subst_sc = {\"g\": 0, \"alpha\": 0, \"mu\": \"mu_sc\", \"V\": 0}\n",
+ " subst_interface = {\"c\": \"c * c_tunnel\", \"alpha\": 0, \"V\": 0}\n",
+ "\n",
+ " templ_sm = discretize(ham, locals=subst_sm, grid_spacing=a)\n",
+ " templ_sc = discretize(ham, locals=subst_sc, grid_spacing=a)\n",
+ " templ_interface = discretize(ham, locals=subst_interface, grid_spacing=a)\n",
+ "\n",
+ " return templ_sm, templ_sc, templ_interface\n",
+ "\n",
+ "\n",
+ "def cylinder_sector(r_out, r_in=0, L=1, L0=0, coverage_angle=360, angle=0, a=10):\n",
+ " \"\"\"Returns the shape function and start coords for a wire with\n",
+ " as cylindrical cross section.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " r_out : int\n",
+ " Outer radius in nm.\n",
+ " r_in : int, optional\n",
+ " Inner radius in nm.\n",
+ " L : int, optional\n",
+ " Length of wire from L0 in nm, -1 if infinite in x-direction.\n",
+ " L0 : int, optional\n",
+ " Start position in x.\n",
+ " coverage_angle : int, optional\n",
+ " Coverage angle in degrees.\n",
+ " angle : int, optional\n",
+ " Angle of tilting from top in degrees.\n",
+ " a : int, optional\n",
+ " Discretization constant in nm.\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " (shape_func, *(start_coords))\n",
+ " \"\"\"\n",
+ " coverage_angle *= np.pi / 360\n",
+ " angle *= np.pi / 180\n",
+ " r_out_sq, r_in_sq = r_out ** 2, r_in ** 2\n",
+ "\n",
+ " def shape(site):\n",
+ " try:\n",
+ " x, y, z = site.pos\n",
+ " except AttributeError:\n",
+ " x, y, z = site\n",
+ " n = (y + 1j * z) * np.exp(1j * angle)\n",
+ " y, z = n.real, n.imag\n",
+ " rsq = y ** 2 + z ** 2\n",
+ " shape_yz = r_in_sq <= rsq < r_out_sq and z >= np.cos(coverage_angle) * np.sqrt(\n",
+ " rsq\n",
+ " )\n",
+ " return (shape_yz and L0 <= x < L) if L > 0 else shape_yz\n",
+ "\n",
+ " r_mid = (r_out + r_in) / 2\n",
+ " start_coords = np.array([L - a, r_mid * np.sin(angle), r_mid * np.cos(angle)])\n",
+ "\n",
+ " return shape, start_coords\n",
+ "\n",
+ "\n",
+ "def is_antisymmetric(H):\n",
+ " return np.allclose(-H, H.T)\n",
+ "\n",
+ "\n",
+ "def cell_mats(lead, params, bias=0):\n",
+ " h = lead.cell_hamiltonian(params=params)\n",
+ " h -= bias * np.identity(len(h))\n",
+ " t = lead.inter_cell_hopping(params=params)\n",
+ " return h, t\n",
+ "\n",
+ "\n",
+ "def get_h_k(lead, params):\n",
+ " h, t = cell_mats(lead, params)\n",
+ "\n",
+ " def h_k(k):\n",
+ " return h + t * np.exp(1j * k) + t.T.conj() * np.exp(-1j * k)\n",
+ "\n",
+ " return h_k\n",
+ "\n",
+ "\n",
+ "def make_skew_symmetric(ham):\n",
+ " \"\"\"\n",
+ " Makes a skew symmetric matrix by a matrix multiplication of a unitary\n",
+ " matrix U. This unitary matrix is taken from the Topology MOOC 0D, but\n",
+ " that is in a different basis. To get to the right basis one multiplies\n",
+ " by [[np.eye(2), 0], [0, sigma_y]].\n",
+ "\n",
+ " Parameters:\n",
+ " -----------\n",
+ " ham : numpy.ndarray\n",
+ " Hamiltonian matrix gotten from sys.cell_hamiltonian()\n",
+ "\n",
+ " Returns:\n",
+ " --------\n",
+ " skew_ham : numpy.ndarray\n",
+ " Skew symmetrized Hamiltonian\n",
+ " \"\"\"\n",
+ " W = ham.shape[0] // 4\n",
+ " I = np.eye(2, dtype=complex)\n",
+ " sigma_y = np.array([[0, 1j], [-1j, 0]], dtype=complex)\n",
+ " U_1 = np.bmat([[I, I], [1j * I, -1j * I]])\n",
+ " U_2 = np.bmat([[I, 0 * I], [0 * I, sigma_y]])\n",
+ " U = U_1 @ U_2\n",
+ " U = np.kron(np.eye(W, dtype=complex), U)\n",
+ " skew_ham = U @ ham @ U.H\n",
+ "\n",
+ " assert is_antisymmetric(skew_ham)\n",
+ "\n",
+ " return skew_ham\n",
+ "\n",
+ "\n",
+ "def calculate_pfaffian(lead, params):\n",
+ " \"\"\"\n",
+ " Calculates the Pfaffian for the infinite system by computing it at k = 0\n",
+ " and k = pi.\n",
+ "\n",
+ " Parameters:\n",
+ " -----------\n",
+ " lead : kwant.builder.InfiniteSystem object\n",
+ " The finalized system.\n",
+ "\n",
+ " \"\"\"\n",
+ " h_k = get_h_k(lead, params)\n",
+ "\n",
+ " skew_h0 = make_skew_symmetric(h_k(0))\n",
+ " skew_h_pi = make_skew_symmetric(h_k(np.pi))\n",
+ "\n",
+ " pf_0 = np.sign(pf.pfaffian(1j * skew_h0, sign_only=True).real)\n",
+ " pf_pi = np.sign(pf.pfaffian(1j * skew_h_pi, sign_only=True).real)\n",
+ " pfaf = pf_0 * pf_pi\n",
+ "\n",
+ " return pfaf\n",
+ "\n",
+ "\n",
+ "def at_interface(site1, site2, shape1, shape2):\n",
+ " return (shape1[0](site1) and shape2[0](site2)) or (\n",
+ " shape2[0](site1) and shape1[0](site2)\n",
+ " )\n",
+ "\n",
+ "\n",
+ "def change_hopping_at_interface(syst, template, shape1, shape2):\n",
+ " for (site1, site2), hop in syst.hopping_value_pairs():\n",
+ " if at_interface(site1, site2, shape1, shape2):\n",
+ " syst[site1, site2] = template[site1, site2]\n",
+ " return syst\n",
+ "\n",
+ "\n",
+ "@memoize\n",
+ "def make_lead(a, r1, r2, coverage_angle, angle, with_shell, which_lead):\n",
+ " \"\"\"Create an infinite cylindrical 3D wire partially covered with a\n",
+ " superconducting (SC) shell.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " a : int\n",
+ " Discretization constant in nm.\n",
+ " r1 : int\n",
+ " Radius of normal part of wire in nm.\n",
+ " r2 : int\n",
+ " Radius of superconductor in nm.\n",
+ " coverage_angle : int\n",
+ " Coverage angle of superconductor in degrees.\n",
+ " angle : int\n",
+ " Angle of tilting of superconductor from top in degrees.\n",
+ " with_shell : bool\n",
+ " Adds shell to the scattering area. If False no SC shell is added and\n",
+ " only a cylindrical wire will be created.\n",
+ " which_lead : str\n",
+ " Name of the potential function of the lead, e.g. `which_lead = 'left'` will\n",
+ " require a function `V_left(z, V_0)` and\n",
+ " `mu_left(mu_func(x, x0, sigma, mu_lead, mu_wire)`.\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " syst : kwant.builder.InfiniteSystem\n",
+ " The finilized kwant system.\n",
+ "\n",
+ " Examples\n",
+ " --------\n",
+ " This doesn't use default parameters because the variables need to be saved,\n",
+ " to a file. So I create a dictionary that is passed to the function.\n",
+ "\n",
+ " >>> syst_params = dict(a=10, angle=0, coverage_angle=185, r1=50,\n",
+ " ... r2=70, with_shell=True)\n",
+ " >>> syst, hopping = make_lead(**syst_params)\n",
+ " \"\"\"\n",
+ "\n",
+ " shape_normal_lead = cylinder_sector(r_out=r1, angle=angle, L=-1, a=a)\n",
+ " shape_sc_lead = cylinder_sector(\n",
+ " r_out=r2, r_in=r1, coverage_angle=coverage_angle, angle=angle, L=-1, a=a\n",
+ " )\n",
+ "\n",
+ " sz = np.array([[1, 0], [0, -1]])\n",
+ " cons_law = np.kron(np.eye(2), -sz)\n",
+ " symmetry = kwant.TranslationalSymmetry((a, 0, 0))\n",
+ " lead = kwant.Builder(\n",
+ " symmetry, conservation_law=cons_law if not with_shell else None\n",
+ " )\n",
+ "\n",
+ " templ_sm, templ_sc, templ_interface = discretized_hamiltonian(\n",
+ " a, which_lead=which_lead\n",
+ " )\n",
+ " templ_sm = apply_peierls_to_template(templ_sm)\n",
+ " lead.fill(templ_sm, *shape_normal_lead)\n",
+ "\n",
+ " if with_shell:\n",
+ " lat = templ_sc.lattice\n",
+ " shape_sc = cylinder_sector(\n",
+ " r_out=r2, r_in=r1, coverage_angle=coverage_angle, angle=angle, L=a, a=a\n",
+ " )\n",
+ "\n",
+ " xyz_offset = get_offset(*shape_sc, lat)\n",
+ "\n",
+ " templ_sc = apply_peierls_to_template(templ_sc, xyz_offset=xyz_offset)\n",
+ " templ_interface = apply_peierls_to_template(templ_interface)\n",
+ " lead.fill(templ_sc, *shape_sc_lead)\n",
+ "\n",
+ " # Adding a tunnel barrier between SM and SC\n",
+ " lead = change_hopping_at_interface(\n",
+ " lead, templ_interface, shape_normal_lead, shape_sc_lead\n",
+ " )\n",
+ "\n",
+ " return lead\n",
+ "\n",
+ "\n",
+ "def apply_peierls_to_template(template, xyz_offset=(0, 0, 0)):\n",
+ " \"\"\"Adds p.orbital argument to the hopping functions.\"\"\"\n",
+ " template = deepcopy(template) # Needed because kwant.Builder is mutable\n",
+ " x0, y0, z0 = xyz_offset\n",
+ " lat = template.lattice\n",
+ " a = np.max(lat.prim_vecs) # lattice contant\n",
+ "\n",
+ " def phase(site1, site2, B_x, B_y, B_z, orbital, e, hbar):\n",
+ " if orbital:\n",
+ " x, y, z = site1.tag\n",
+ " direction = site1.tag - site2.tag\n",
+ " A = [B_y * (z - z0) - B_z * (y - y0), 0, B_x * (y - y0)]\n",
+ " A = np.dot(A, direction) * a ** 2 * 1e-18 * e / hbar\n",
+ " phase = np.exp(-1j * A)\n",
+ " if lat.norbs == 2: # No PH degrees of freedom\n",
+ " return phase\n",
+ " elif lat.norbs == 4:\n",
+ " return np.array(\n",
+ " [phase, phase.conj(), phase, phase.conj()], dtype=\"complex128\"\n",
+ " )\n",
+ " else: # No orbital phase\n",
+ " return 1\n",
+ "\n",
+ " for (site1, site2), hop in template.hopping_value_pairs():\n",
+ " template[site1, site2] = combine(hop, phase, operator.mul, 2)\n",
+ " return template\n",
+ "\n",
+ "\n",
+ "def get_offset(shape, start, lat):\n",
+ " a = np.max(lat.prim_vecs)\n",
+ " coords = [site.pos for site in lat.shape(shape, start)()]\n",
+ " xyz_offset = np.mean(coords, axis=0)\n",
+ " return xyz_offset\n",
+ "\n",
+ "\n",
+ "def translation_ev(h, t, tol=1e6):\n",
+ " \"\"\"Compute the eigenvalues of the translation operator of a lead.\n",
+ "\n",
+ " Adapted from kwant.physics.leads.modes.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " h : numpy array, real or complex, shape (N, N) The unit cell\n",
+ " Hamiltonian of the lead unit cell.\n",
+ " t : numpy array, real or complex, shape (N, M)\n",
+ " The hopping matrix from a lead cell to the one on which self-energy\n",
+ " has to be calculated (and any other hopping in the same direction).\n",
+ " tol : float\n",
+ " Numbers and differences are considered zero when they are smaller\n",
+ " than `tol` times the machine precision.\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " ev : numpy array\n",
+ " Eigenvalues of the translation operator in the form lambda=r*exp(i*k),\n",
+ " for |r|=1 they are propagating modes.\n",
+ " \"\"\"\n",
+ " a, b = kwant.physics.leads.setup_linsys(h, t, tol, None).eigenproblem\n",
+ " ev = kwant.physics.leads.unified_eigenproblem(a, b, tol=tol)[0]\n",
+ " return ev\n",
+ "\n",
+ "\n",
+ "def gap_minimizer(lead, params, energy):\n",
+ " \"\"\"Function that minimizes a function to find the band gap.\n",
+ " This objective function checks if there are progagating modes at a\n",
+ " certain energy. Returns zero if there is a propagating mode.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " lead : kwant.builder.InfiniteSystem object\n",
+ " The finalized infinite system.\n",
+ " params : dict\n",
+ " A dict that is used to store Hamiltonian parameters.\n",
+ " energy : float\n",
+ " Energy at which this function checks for propagating modes.\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " minimized_scalar : float\n",
+ " Value that is zero when there is a propagating mode.\n",
+ " \"\"\"\n",
+ " h, t = cell_mats(lead, params, bias=energy)\n",
+ " ev = translation_ev(h, t)\n",
+ " norm = (ev * ev.conj()).real\n",
+ " return np.min(np.abs(norm - 1))\n",
+ "\n",
+ "\n",
+ "def gap_from_modes(lead, params, tol=1e-6):\n",
+ " \"\"\"Finds the gapsize by peforming a binary search of the modes with a\n",
+ " tolarance of tol.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " lead : kwant.builder.InfiniteSystem object\n",
+ " The finalized infinite system.\n",
+ " params : dict\n",
+ " A dict that is used to store Hamiltonian parameters.\n",
+ " tol : float\n",
+ " The precision of the binary search.\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " gap : float\n",
+ " Size of the gap.\n",
+ "\n",
+ " Notes\n",
+ " -----\n",
+ " For use with `lead = funcs.make_lead()`.\n",
+ " \"\"\"\n",
+ " Es = kwant.physics.Bands(lead, params=params)(k=0)\n",
+ " lim = [0, np.abs(Es).min()]\n",
+ " if gap_minimizer(lead, params, energy=0) < 1e-15:\n",
+ " # No band gap\n",
+ " gap = 0\n",
+ " else:\n",
+ " while lim[1] - lim[0] > tol:\n",
+ " energy = sum(lim) / 2\n",
+ " par = gap_minimizer(lead, params, energy)\n",
+ " if par < 1e-10:\n",
+ " lim[1] = energy\n",
+ " else:\n",
+ " lim[0] = energy\n",
+ " gap = sum(lim) / 2\n",
+ " return gap"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import adaptive\n",
+ "adaptive.notebook_extension()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "\n",
+ "import phase_diagram as funcs\n",
+ "\n",
+ "lead_pars = dict(\n",
+ " a=10, r1=50, r2=70, coverage_angle=135, angle=45, with_shell=True, which_lead=''\n",
+ ")\n",
+ "\n",
+ "params = dict(\n",
+ " alpha=20,\n",
+ " B_x=0,\n",
+ " B_y=0,\n",
+ " B_z=0,\n",
+ " Delta=110,\n",
+ " g=50,\n",
+ " orbital=True,\n",
+ " mu_sc=100,\n",
+ " c_tunnel=3 / 4,\n",
+ " V_r=-50,\n",
+ " intrinsic_sc=False,\n",
+ " mu_=lambda x0, sigma, mu_lead, mu_wire: mu_lead,\n",
+ " V_=lambda z, V_0, V_r, V_l, x0, sigma, r1: 0,\n",
+ " V_0=None,\n",
+ " V_l=None,\n",
+ " mu_lead=10,\n",
+ " mu_wire=None,\n",
+ " r1=None,\n",
+ " sigma=None,\n",
+ " x0=None,\n",
+ " **funcs.constants.__dict__\n",
+ ")\n",
+ "\n",
+ "def pf(xy, params=params, lead_pars=lead_pars):\n",
+ " import phase_diagram as funcs\n",
+ " params[\"B_x\"], params[\"mu_lead\"] = xy\n",
+ " lead = funcs.make_lead(**lead_pars).finalized()\n",
+ " return funcs.calculate_pfaffian(lead, params)\n",
+ "\n",
+ "def smallest_gap(xy, params=params, lead_pars=lead_pars):\n",
+ " import phase_diagram as funcs\n",
+ " params[\"B_x\"], params[\"mu_lead\"] = xy\n",
+ " lead = funcs.make_lead(**lead_pars).finalized()\n",
+ " pf = funcs.calculate_pfaffian(lead, params)\n",
+ " gap = funcs.gap_from_modes(lead, params)\n",
+ " return pf * gap "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "learner = adaptive.Learner2D(smallest_gap, bounds=[(0, 2), (0, 35)])\n",
+ "# learner.load('phase_diagram.pickle')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "runner = adaptive.Runner(learner, goal=lambda l: l.npoints > 60000)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "runner.live_info()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%%output size=100\n",
+ "learner.plot(tri_alpha=0.4, n=None)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "learner.save('phase_diagram.pickle')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "language_info": {
+ "name": "python",
+ "pygments_lexer": "ipython3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/phase_diagram.py b/phase_diagram.py
new file mode 100644
index 0000000000000000000000000000000000000000..f154fca12acc4b54ccd89a21fcbd484276e76804
--- /dev/null
+++ b/phase_diagram.py
@@ -0,0 +1,472 @@
+
+
+import inspect
+import operator
+import sys
+from collections import OrderedDict
+from copy import copy, deepcopy
+from functools import wraps
+from itertools import product
+from types import SimpleNamespace
+
+import kwant
+import numpy as np
+import scipy.constants
+from kwant.continuum.discretizer import discretize
+from skimage import measure
+
+import adaptive
+import pfaffian as pf
+
+assert sys.version_info >= (3, 6), "Use Python ≥3.6"
+
+# Parameters taken from arXiv:1204.2792
+# All constant parameters, mostly fundamental constants, in a SimpleNamespace.
+constants = SimpleNamespace(
+ m_eff=0.015 * scipy.constants.m_e, # effective mass in kg
+ hbar=scipy.constants.hbar,
+ m_e=scipy.constants.m_e,
+ eV=scipy.constants.eV,
+ e=scipy.constants.e,
+ c=1e18 / (scipy.constants.eV * 1e-3), # to get to meV * nm^2
+ mu_B=scipy.constants.physical_constants["Bohr magneton in eV/T"][0] * 1e3,
+)
+
+constants.t = (constants.hbar ** 2 / (2 * constants.m_eff)) * constants.c
+
+
+def get_names(sig):
+ names = [
+ (name, value)
+ for name, value in sig.parameters.items()
+ if value.kind
+ in (inspect.Parameter.POSITIONAL_OR_KEYWORD, inspect.Parameter.KEYWORD_ONLY)
+ ]
+ return OrderedDict(names)
+
+
+def filter_kwargs(sig, names, kwargs):
+ names_in_kwargs = [(name, value) for name, value in kwargs.items() if name in names]
+ return OrderedDict(names_in_kwargs)
+
+
+def skip_pars(names1, names2, num_skipped):
+ skipped_pars1 = list(names1.keys())[:num_skipped]
+ skipped_pars2 = list(names2.keys())[:num_skipped]
+ if skipped_pars1 == skipped_pars2:
+ pars1 = list(names1.values())[num_skipped:]
+ pars2 = list(names2.values())[num_skipped:]
+ else:
+ raise Exception("First {} arguments " "have to be the same".format(num_skipped))
+ return pars1, pars2
+
+
+def combine(f, g, operator, num_skipped=0):
+ if not callable(f) or not callable(g):
+ raise Exception("One of the functions is not a function")
+
+ sig1 = inspect.signature(f)
+ sig2 = inspect.signature(g)
+
+ names1 = get_names(sig1)
+ names2 = get_names(sig2)
+
+ pars1, pars2 = skip_pars(names1, names2, num_skipped)
+ skipped_pars = list(names1.values())[:num_skipped]
+
+ pars1_names = {p.name for p in pars1}
+ pars2 = [p for p in pars2 if p.name not in pars1_names]
+
+ parameters = pars1 + pars2
+ kind = inspect.Parameter.POSITIONAL_OR_KEYWORD
+ parameters = [p.replace(kind=kind) for p in parameters]
+ parameters = skipped_pars + parameters
+
+ def wrapped(*args):
+ d = {p.name: arg for arg, p in zip(args, parameters)}
+ fval = f(*[d[name] for name in names1.keys()])
+ gval = g(*[d[name] for name in names2.keys()])
+ return operator(fval, gval)
+
+ wrapped.__signature__ = inspect.Signature(parameters=parameters)
+ return wrapped
+
+
+def memoize(obj):
+ cache = obj.cache = {}
+
+ @wraps(obj)
+ def memoizer(*args, **kwargs):
+ key = str(args) + str(kwargs)
+ if key not in cache:
+ cache[key] = obj(*args, **kwargs)
+ return cache[key]
+
+ return memoizer
+
+
+@memoize
+def discretized_hamiltonian(a, which_lead=None):
+ ham = (
+ "(0.5 * hbar**2 * (k_x**2 + k_y**2 + k_z**2) / m_eff * c - mu + V) * kron(sigma_0, sigma_z) + "
+ "alpha * (k_y * kron(sigma_x, sigma_z) - k_x * kron(sigma_y, sigma_z)) + "
+ "0.5 * g * mu_B * (B_x * kron(sigma_x, sigma_0) + B_y * kron(sigma_y, sigma_0) + B_z * kron(sigma_z, sigma_0)) + "
+ "Delta * kron(sigma_0, sigma_x)"
+ )
+
+ subst_sm = {"Delta": 0}
+
+ if which_lead is not None:
+ subst_sm["V"] = f"V_{which_lead}(z, V_0, V_r, V_l, x0, sigma, r1)"
+ subst_sm["mu"] = f"mu_{which_lead}(x0, sigma, mu_lead, mu_wire)"
+ else:
+ subst_sm["V"] = "V(x, z, V_0, V_r, V_l, x0, sigma, r1)"
+ subst_sm["mu"] = "mu(x, x0, sigma, mu_lead, mu_wire)"
+
+ subst_sc = {"g": 0, "alpha": 0, "mu": "mu_sc", "V": 0}
+ subst_interface = {"c": "c * c_tunnel", "alpha": 0, "V": 0}
+
+ templ_sm = discretize(ham, locals=subst_sm, grid_spacing=a)
+ templ_sc = discretize(ham, locals=subst_sc, grid_spacing=a)
+ templ_interface = discretize(ham, locals=subst_interface, grid_spacing=a)
+
+ return templ_sm, templ_sc, templ_interface
+
+
+def cylinder_sector(r_out, r_in=0, L=1, L0=0, coverage_angle=360, angle=0, a=10):
+ """Returns the shape function and start coords for a wire with
+ as cylindrical cross section.
+
+ Parameters
+ ----------
+ r_out : int
+ Outer radius in nm.
+ r_in : int, optional
+ Inner radius in nm.
+ L : int, optional
+ Length of wire from L0 in nm, -1 if infinite in x-direction.
+ L0 : int, optional
+ Start position in x.
+ coverage_angle : int, optional
+ Coverage angle in degrees.
+ angle : int, optional
+ Angle of tilting from top in degrees.
+ a : int, optional
+ Discretization constant in nm.
+
+ Returns
+ -------
+ (shape_func, *(start_coords))
+ """
+ coverage_angle *= np.pi / 360
+ angle *= np.pi / 180
+ r_out_sq, r_in_sq = r_out ** 2, r_in ** 2
+
+ def shape(site):
+ try:
+ x, y, z = site.pos
+ except AttributeError:
+ x, y, z = site
+ n = (y + 1j * z) * np.exp(1j * angle)
+ y, z = n.real, n.imag
+ rsq = y ** 2 + z ** 2
+ shape_yz = r_in_sq <= rsq < r_out_sq and z >= np.cos(coverage_angle) * np.sqrt(
+ rsq
+ )
+ return (shape_yz and L0 <= x < L) if L > 0 else shape_yz
+
+ r_mid = (r_out + r_in) / 2
+ start_coords = np.array([L - a, r_mid * np.sin(angle), r_mid * np.cos(angle)])
+
+ return shape, start_coords
+
+
+def is_antisymmetric(H):
+ return np.allclose(-H, H.T)
+
+
+def cell_mats(lead, params, bias=0):
+ h = lead.cell_hamiltonian(params=params)
+ h -= bias * np.identity(len(h))
+ t = lead.inter_cell_hopping(params=params)
+ return h, t
+
+
+def get_h_k(lead, params):
+ h, t = cell_mats(lead, params)
+
+ def h_k(k):
+ return h + t * np.exp(1j * k) + t.T.conj() * np.exp(-1j * k)
+
+ return h_k
+
+
+def make_skew_symmetric(ham):
+ """
+ Makes a skew symmetric matrix by a matrix multiplication of a unitary
+ matrix U. This unitary matrix is taken from the Topology MOOC 0D, but
+ that is in a different basis. To get to the right basis one multiplies
+ by [[np.eye(2), 0], [0, sigma_y]].
+
+ Parameters:
+ -----------
+ ham : numpy.ndarray
+ Hamiltonian matrix gotten from sys.cell_hamiltonian()
+
+ Returns:
+ --------
+ skew_ham : numpy.ndarray
+ Skew symmetrized Hamiltonian
+ """
+ W = ham.shape[0] // 4
+ I = np.eye(2, dtype=complex)
+ sigma_y = np.array([[0, 1j], [-1j, 0]], dtype=complex)
+ U_1 = np.bmat([[I, I], [1j * I, -1j * I]])
+ U_2 = np.bmat([[I, 0 * I], [0 * I, sigma_y]])
+ U = U_1 @ U_2
+ U = np.kron(np.eye(W, dtype=complex), U)
+ skew_ham = U @ ham @ U.H
+
+ assert is_antisymmetric(skew_ham)
+
+ return skew_ham
+
+
+def calculate_pfaffian(lead, params):
+ """
+ Calculates the Pfaffian for the infinite system by computing it at k = 0
+ and k = pi.
+
+ Parameters:
+ -----------
+ lead : kwant.builder.InfiniteSystem object
+ The finalized system.
+
+ """
+ h_k = get_h_k(lead, params)
+
+ skew_h0 = make_skew_symmetric(h_k(0))
+ skew_h_pi = make_skew_symmetric(h_k(np.pi))
+
+ pf_0 = np.sign(pf.pfaffian(1j * skew_h0, sign_only=True).real)
+ pf_pi = np.sign(pf.pfaffian(1j * skew_h_pi, sign_only=True).real)
+ pfaf = pf_0 * pf_pi
+
+ return pfaf
+
+
+def at_interface(site1, site2, shape1, shape2):
+ return (shape1[0](site1) and shape2[0](site2)) or (
+ shape2[0](site1) and shape1[0](site2)
+ )
+
+
+def change_hopping_at_interface(syst, template, shape1, shape2):
+ for (site1, site2), hop in syst.hopping_value_pairs():
+ if at_interface(site1, site2, shape1, shape2):
+ syst[site1, site2] = template[site1, site2]
+ return syst
+
+
+@memoize
+def make_lead(a, r1, r2, coverage_angle, angle, with_shell, which_lead):
+ """Create an infinite cylindrical 3D wire partially covered with a
+ superconducting (SC) shell.
+
+ Parameters
+ ----------
+ a : int
+ Discretization constant in nm.
+ r1 : int
+ Radius of normal part of wire in nm.
+ r2 : int
+ Radius of superconductor in nm.
+ coverage_angle : int
+ Coverage angle of superconductor in degrees.
+ angle : int
+ Angle of tilting of superconductor from top in degrees.
+ with_shell : bool
+ Adds shell to the scattering area. If False no SC shell is added and
+ only a cylindrical wire will be created.
+ which_lead : str
+ Name of the potential function of the lead, e.g. `which_lead = 'left'` will
+ require a function `V_left(z, V_0)` and
+ `mu_left(mu_func(x, x0, sigma, mu_lead, mu_wire)`.
+
+ Returns
+ -------
+ syst : kwant.builder.InfiniteSystem
+ The finilized kwant system.
+
+ Examples
+ --------
+ This doesn't use default parameters because the variables need to be saved,
+ to a file. So I create a dictionary that is passed to the function.
+
+ >>> syst_params = dict(a=10, angle=0, coverage_angle=185, r1=50,
+ ... r2=70, with_shell=True)
+ >>> syst, hopping = make_lead(**syst_params)
+ """
+
+ shape_normal_lead = cylinder_sector(r_out=r1, angle=angle, L=-1, a=a)
+ shape_sc_lead = cylinder_sector(
+ r_out=r2, r_in=r1, coverage_angle=coverage_angle, angle=angle, L=-1, a=a
+ )
+
+ sz = np.array([[1, 0], [0, -1]])
+ cons_law = np.kron(np.eye(2), -sz)
+ symmetry = kwant.TranslationalSymmetry((a, 0, 0))
+ lead = kwant.Builder(
+ symmetry, conservation_law=cons_law if not with_shell else None
+ )
+
+ templ_sm, templ_sc, templ_interface = discretized_hamiltonian(
+ a, which_lead=which_lead
+ )
+ templ_sm = apply_peierls_to_template(templ_sm)
+ lead.fill(templ_sm, *shape_normal_lead)
+
+ if with_shell:
+ lat = templ_sc.lattice
+ shape_sc = cylinder_sector(
+ r_out=r2, r_in=r1, coverage_angle=coverage_angle, angle=angle, L=a, a=a
+ )
+
+ xyz_offset = get_offset(*shape_sc, lat)
+
+ templ_sc = apply_peierls_to_template(templ_sc, xyz_offset=xyz_offset)
+ templ_interface = apply_peierls_to_template(templ_interface)
+ lead.fill(templ_sc, *shape_sc_lead)
+
+ # Adding a tunnel barrier between SM and SC
+ lead = change_hopping_at_interface(
+ lead, templ_interface, shape_normal_lead, shape_sc_lead
+ )
+
+ return lead
+
+
+def apply_peierls_to_template(template, xyz_offset=(0, 0, 0)):
+ """Adds p.orbital argument to the hopping functions."""
+ template = deepcopy(template) # Needed because kwant.Builder is mutable
+ x0, y0, z0 = xyz_offset
+ lat = template.lattice
+ a = np.max(lat.prim_vecs) # lattice contant
+
+ def phase(site1, site2, B_x, B_y, B_z, orbital, e, hbar):
+ if orbital:
+ x, y, z = site1.tag
+ direction = site1.tag - site2.tag
+ A = [B_y * (z - z0) - B_z * (y - y0), 0, B_x * (y - y0)]
+ A = np.dot(A, direction) * a ** 2 * 1e-18 * e / hbar
+ phase = np.exp(-1j * A)
+ if lat.norbs == 2: # No PH degrees of freedom
+ return phase
+ elif lat.norbs == 4:
+ return np.array(
+ [phase, phase.conj(), phase, phase.conj()], dtype="complex128"
+ )
+ else: # No orbital phase
+ return 1
+
+ for (site1, site2), hop in template.hopping_value_pairs():
+ template[site1, site2] = combine(hop, phase, operator.mul, 2)
+ return template
+
+
+def get_offset(shape, start, lat):
+ a = np.max(lat.prim_vecs)
+ coords = [site.pos for site in lat.shape(shape, start)()]
+ xyz_offset = np.mean(coords, axis=0)
+ return xyz_offset
+
+
+def translation_ev(h, t, tol=1e6):
+ """Compute the eigenvalues of the translation operator of a lead.
+
+ Adapted from kwant.physics.leads.modes.
+
+ Parameters
+ ----------
+ h : numpy array, real or complex, shape (N, N) The unit cell
+ Hamiltonian of the lead unit cell.
+ t : numpy array, real or complex, shape (N, M)
+ The hopping matrix from a lead cell to the one on which self-energy
+ has to be calculated (and any other hopping in the same direction).
+ tol : float
+ Numbers and differences are considered zero when they are smaller
+ than `tol` times the machine precision.
+
+ Returns
+ -------
+ ev : numpy array
+ Eigenvalues of the translation operator in the form lambda=r*exp(i*k),
+ for |r|=1 they are propagating modes.
+ """
+ a, b = kwant.physics.leads.setup_linsys(h, t, tol, None).eigenproblem
+ ev = kwant.physics.leads.unified_eigenproblem(a, b, tol=tol)[0]
+ return ev
+
+
+def gap_minimizer(lead, params, energy):
+ """Function that minimizes a function to find the band gap.
+ This objective function checks if there are progagating modes at a
+ certain energy. Returns zero if there is a propagating mode.
+
+ Parameters
+ ----------
+ lead : kwant.builder.InfiniteSystem object
+ The finalized infinite system.
+ params : dict
+ A dict that is used to store Hamiltonian parameters.
+ energy : float
+ Energy at which this function checks for propagating modes.
+
+ Returns
+ -------
+ minimized_scalar : float
+ Value that is zero when there is a propagating mode.
+ """
+ h, t = cell_mats(lead, params, bias=energy)
+ ev = translation_ev(h, t)
+ norm = (ev * ev.conj()).real
+ return np.min(np.abs(norm - 1))
+
+
+def gap_from_modes(lead, params, tol=1e-6):
+ """Finds the gapsize by peforming a binary search of the modes with a
+ tolarance of tol.
+
+ Parameters
+ ----------
+ lead : kwant.builder.InfiniteSystem object
+ The finalized infinite system.
+ params : dict
+ A dict that is used to store Hamiltonian parameters.
+ tol : float
+ The precision of the binary search.
+
+ Returns
+ -------
+ gap : float
+ Size of the gap.
+
+ Notes
+ -----
+ For use with `lead = funcs.make_lead()`.
+ """
+ Es = kwant.physics.Bands(lead, params=params)(k=0)
+ lim = [0, np.abs(Es).min()]
+ if gap_minimizer(lead, params, energy=0) < 1e-15:
+ # No band gap
+ gap = 0
+ else:
+ while lim[1] - lim[0] > tol:
+ energy = sum(lim) / 2
+ par = gap_minimizer(lead, params, energy)
+ if par < 1e-10:
+ lim[1] = energy
+ else:
+ lim[0] = energy
+ gap = sum(lim) / 2
+ return gap