From 33980a70d31a24ead9ebde2ae87f673a5505f8f5 Mon Sep 17 00:00:00 2001
From: Bas Nijholt <basnijholt@gmail.com>
Date: Mon, 23 Sep 2019 13:28:43 +0200
Subject: [PATCH] put code in module
---
phase_diagram.ipynb | 497 +-------------------------------------------
1 file changed, 6 insertions(+), 491 deletions(-)
diff --git a/phase_diagram.ipynb b/phase_diagram.ipynb
index 4b73362..009d5c9 100644
--- a/phase_diagram.ipynb
+++ b/phase_diagram.ipynb
@@ -1,486 +1,5 @@
{
"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,
@@ -497,11 +16,10 @@
"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",
+ " a=10, r1=50, r2=70, coverage_angle=135, angle=45, with_shell=True, which_lead=\"\"\n",
")\n",
"\n",
"params = dict(\n",
@@ -528,19 +46,23 @@
" **funcs.constants.__dict__\n",
")\n",
"\n",
+ "\n",
"def pf(xy, params=params, lead_pars=lead_pars):\n",
" import phase_diagram as funcs\n",
+ "\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",
+ "\n",
"def smallest_gap(xy, params=params, lead_pars=lead_pars):\n",
" import phase_diagram as funcs\n",
+ "\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 "
+ " return pf * gap"
]
},
{
@@ -571,13 +93,6 @@
"runner.live_info()"
]
},
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- },
{
"cell_type": "code",
"execution_count": null,
--
GitLab