In [1]:
import sys

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

init_notebook()

import scipy
from time import time
import scipy.linalg as sla

randn = np.random.randn


def chern_torus(W=5, L=5):

    disorder_potential = randn(W, L)

    def shape(pos):
        (x, y) = pos
        return (0 <= x < W) * (0 <= y < L)

    def onsite(site, p):
        t, mu, dis_amplitude = p.t, p.mu, p.disorder
        (x, y) = site.pos
        return (4 * t - mu) * pauli.sz

    def hopx(site1, site2, p):
        t, delta = p.t, p.delta
        return -t * pauli.sz + 1j * p.delta * pauli.sx

    def hopy(site1, site2, p):
        t, delta = p.t, p.delta
        return -t * pauli.sz - 1j * p.delta * pauli.sy

    lat = kwant.lattice.square()
    syst = kwant.Builder()
    syst[lat.shape(shape, (0, 0))] = onsite
    syst[kwant.HoppingKind((1, 0), lat)] = hopx
    syst[kwant.HoppingKind((0, 1), lat)] = hopy
    syst[kwant.HoppingKind((-W + 1, 0), lat)] = hopx
    syst[kwant.HoppingKind((0, -L + 1), lat)] = hopy
    return syst.finalized()


def projected_operators(syst, p, energy=0):
    x, y = np.array([i.pos for i in syst.sites]).T
    # Double all entries to take orbitals into account.
    x = np.resize(x, (2, len(x))).T.flatten()
    y = np.resize(y, (2, len(y))).T.flatten()
    x *= 2 * np.pi / (np.max(x) + 1)
    y *= 2 * np.pi / (np.max(y) + 1)
    op_x = np.diag(np.exp(1j * x))
    op_y = np.diag(np.exp(1j * y))

    ham = syst.hamiltonian_submatrix(params=dict(p=p))
    ham -= 0 * np.identity(len(ham))
    energies, states = np.linalg.eigh(ham)

    projector = states[:, energies < energy]

    op_x = projector.T.conj() @ op_x @ projector
    op_y = projector.T.conj() @ op_y @ projector
    return op_x, op_y


def two_terminal(L, W):
    t = 1.0

    def shape(pos):
        (x, y) = pos
        return 0 <= y < W and 0 <= x < L

    def lead_shape(pos):
        (x, y) = pos
        return 0 <= y < W

    lat = kwant.lattice.square()
    syst = kwant.Builder()

    # definition of system
    syst[lat.shape(shape, (0, 0))] = 3 * t
    syst[kwant.HoppingKind((1, 0), lat)] = -t
    syst[kwant.HoppingKind((0, 1), lat)] = -t

    # definition of leads
    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
    lead[lat.shape(lead_shape, (0, 0))] = 3 * t
    lead[kwant.HoppingKind((1, 0), lat)] = -t
    lead[kwant.HoppingKind((0, 1), lat)] = -t

    # attaching leads
    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())

    return syst.finalized()


def diag_time(N):
    syst = two_terminal(N, N)

    start = time()
    ham = syst.hamiltonian_submatrix()
    ev, evec = sla.eigh(ham)
    res = time() - start
    return res


def smat_time(N):
    syst = two_terminal(N, N)

    start = time()

    smat = kwant.smatrix(syst)

    res = time() - start
    return res
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:04.747821.