In [1]:
import sys

sys.path.append("../../code")
from init_mooc_nb import *

init_notebook()
%output size = 150
import scipy
from matplotlib import cm

bhz_parameters = {
    "topo": {"A": 0.5, "B": 1.00, "D": 0.0, "M": 1.0, "del_z": 0.0},
    "triv": {"A": 0.5, "B": 1.00, "D": 0.0, "M": -1.0, "del_z": 0.0},
    "topo2": {"A": 0.5, "B": 1.00, "D": 0.3, "M": 1.0, "del_z": 0.0},
    "slowed": {"A": 0.05, "B": 0.08, "D": 0.15, "M": -0.3, "del_z": 0.5},
}

# Onsite and hoppings for bhz model


def onsite(site, p):
    return (p.M - 4 * p.B) * pauli.s0sz - 4 * p.D * pauli.s0s0 + p.del_z * pauli.sysy


def hopx(site1, site2, p):
    return p.B * pauli.s0sz + p.D * pauli.s0s0 + 1j * p.A * pauli.szsx


def hopy(site1, site2, p):
    return p.B * pauli.s0sz + p.D * pauli.s0s0 - 1j * p.A * pauli.s0sy


def bhz(w=20):
    """ Make ribbon system with bhz model.

    slowed parameters are used on the edge for finite size system.
    """
    lat = kwant.lattice.square()

    def hopping_x(site1, site2, p):
        x1, y1 = site1.pos
        x2, y2 = site2.pos
        return hopx(site1, site2, p) * np.exp(-0.5j * p.Bz * (x1 - x2) * (y1 + y2))

    slowed_par = SimpleNamespace(Bz=0, **bhz_parameters["slowed"])
    if w is None:
        syst = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))
        syst[lat.shape(lambda pos: True, (0, 0))] = onsite
        syst[kwant.HoppingKind((1, 0), lat)] = hopping_x
        syst[kwant.HoppingKind((0, 1), lat)] = hopy
    else:
        syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
        syst[(lat(0, i) for i in range(w))] = onsite
        syst[kwant.HoppingKind((1, 0), lat)] = hopping_x
        syst[kwant.HoppingKind((0, 1), lat)] = hopy

        syst[lat(0, -1)] = lambda site, p: onsite(site, slowed_par)
        syst[lat(1, -1), lat(0, -1)] = lambda site1, site2, p: hopping_x(
            site1, site2, slowed_par
        )
        syst[lat(0, 0), lat(0, -1)] = hopy

        syst[lat(0, w)] = lambda site, p: onsite(site, slowed_par)
        syst[lat(1, w), lat(0, w)] = lambda site1, site2, p: hopping_x(
            site1, site2, slowed_par
        )
        syst[lat(0, w), lat(0, w - 1)] = hopy

    return syst


def bhz_cylinder(w=3):
    """ Make cylinder system with bhz model. """

    def ribbon_shape(pos):
        (x, y) = pos
        return 0 <= y < w

    lat = kwant.lattice.square(norbs=4)
    sym = kwant.TranslationalSymmetry((1, 0))
    syst = kwant.Builder(sym)

    def hopping_x(site1, site2, p):
        x1, y1 = site1.pos
        x2, y2 = site2.pos
        return hopx(site1, site2, p) * np.exp(-0.5j * p.Bz * (x1 - x2) * (y1 + y2))

    def hopy_phase(site1, site2, p):
        x1, y1 = site1.pos
        x2, y2 = site2.pos
        return hopy(site1, site2, p) * np.exp(1j * p.ky)

    syst[lat.shape(ribbon_shape, (0, 0))] = onsite
    syst[kwant.HoppingKind((1, 0), lat)] = hopping_x
    syst[kwant.HoppingKind((0, 1), lat)] = hopy
    syst[kwant.HoppingKind((0, -w + 1), lat)] = hopy_phase

    syst[lat(0, 0)] = onsite
    syst[lat(0, w - 1)] = onsite

    return syst


def make_lead(t, trs=None):
    def ribbon_shape(pos):
        (x, y) = pos
        return 0 <= y < 2

    lat = kwant.lattice.square(norbs=4)
    sym = kwant.TranslationalSymmetry((1, 0))
    syst = kwant.Builder(sym, time_reversal=1j * pauli.sys0)

    syst[lat.shape(ribbon_shape, (0, 0))] = 1.8 * t * pauli.s0sz
    syst[kwant.HoppingKind((1, 0), lat)] = -t * pauli.s0sz
    return syst


def make_scatter_sys():
    def shape(pos):
        x, y = pos
        return (x == 0) * (0 <= y < 3)

    lat = kwant.lattice.square(norbs=4)
    syst = kwant.Builder()
    syst[lat.shape(shape, (0, 0))] = onsite
    syst[kwant.HoppingKind((0, 1), lat)] = hopy

    lead_cylinder = bhz_cylinder()
    lead = make_lead(1.0)
    syst.attach_lead(lead.reversed())
    syst.attach_lead(lead_cylinder)
    syst = syst.finalized()
    return syst


def scattering_det_pfaff(syst, p):
    def pfaffian(syst, p, ky):
        p.ky = ky
        smat = kwant.smatrix(syst, energy=0.0, params=dict(p=p)).data
        # since we get relatively large numerical errors we project the matrix on
        # the space of antisymmetric matrices
        smat = 0.5 * (smat - smat.T)
        return pf.pfaffian(smat)

    pfaff = [pfaffian(syst, p, 0), pfaffian(syst, p, np.pi)]

    ks = np.linspace(0.0, np.pi, 50)
    det = [np.linalg.det(kwant.smatrix(syst, energy=0.0, params=dict(p=p)).data) for p.ky in ks]
    det = np.array(det)

    phase = np.angle(pfaff[0]) + 0.5 * np.cumsum(np.angle(det[1:] / det[:-1]))
    kdims = ["$k_y$", "phase"]
    plot = holoviews.Path((ks[1:], phase), kdims=kdims).opts(style={"color": "b"})
    plot *= holoviews.Points(([0, np.pi], np.angle(pfaff)), kdims=kdims).opts(
        style={"color": "g"}
    )
    xlims, ylims = slice(-0.2, np.pi + 0.2), slice(-np.pi - 0.2, np.pi + 0.2)
    pi_ticks = [(-np.pi, r"$-\pi$"), (0, "$0$"), (np.pi, r"$\pi$")]
    ticks = {"xticks": [(0, "0"), (np.pi, "$\pi$")], "yticks": pi_ticks}
    return plot.relabel("Winding", depth=1)[xlims, ylims].opts(plot=ticks)


def title(p):
    title = r"$A={:.2}$, $B={:.2}$, $D={:.2}$, $M={:.2}$"
    return title.format(p.A, p.B, p.D, p.M)
Populated the namespace with:
np, matplotlib, kwant, holoviews, init_notebook, SimpleNamespace, pprint_matrix, scientific_number, pretty_fmt_complex, plt, pf, display_html
from code/edx_components:
MoocVideo, MoocDiscussion, MoocCheckboxesAssessment, MoocMultipleChoiceAssessment, MoocSelfAssessment
from code/functions:
spectrum, hamiltonian_array, h_k, pauli
Using kwant 1.4.2 and holoviews 1.13.2
Executed on 2020-06-07 at 08:22:31.860422.