In [2]:
%run ../code/init_mooc_nb.ipy

import scipy
from matplotlib import cm
from matplotlib import gridspec  
import tinyarray as ta
import scipy.sparse.linalg as sl


sigma0 = np.array([[1, 0], [0, 1]])
sigmax = np.array([[0, 1], [1, 0]])
sigmay = np.array([[0, -1j], [1j, 0]])
sigmaz = np.array([[1, 0], [0, -1]])


class SimpleNamespace(object):
    def __init__(self, **kwargs):
        self.__dict__.update(kwargs)

        
# Onsite and hoppings matrices used for building BHZ model
def onsite(site, par):
    A1, A2, B1, B2, D1, D2, C, M = par.A1, par.A2, par.B1, par.B2, par.D1, par.D2, par.C, par.M
    return (C + 2 * D1 + 4 * D2) * np.kron(sigma0, sigma0) + (M + 2 * B1 + 4 * B2) * np.kron(sigma0, sigmaz)


def hopx(site1, site2, par):
    A1, A2, B1, B2, D1, D2, C, M = par.A1, par.A2, par.B1, par.B2, par.D1, par.D2, par.C, par.M
    return - D2 * np.kron(sigma0, sigma0) - B2 * np.kron(sigma0, sigmaz) + A2 * 0.5j * np.kron(sigmax, sigmax)


def hopy(site1, site2, par):
    A1, A2, B1, B2, D1, D2, C, M = par.A1, par.A2, par.B1, par.B2, par.D1, par.D2, par.C, par.M
    return - D2 * np.kron(sigma0, sigma0) - B2 * np.kron(sigma0, sigmaz) + A2 * 0.5j * np.kron(sigmay, sigmax)


def hopz(site1, site2, par):
    A1, A2, B1, B2, D1, D2, C, M = par.A1, par.A2, par.B1, par.B2, par.D1, par.D2, par.C, par.M
    return - D1 * np.kron(sigma0, sigma0) - B1 * np.kron(sigma0, sigmaz) + A1 * 0.5j * np.kron(sigmaz, sigmax)

    
def bhz_slab(L, W, H):
    #lat = kwant.lattice.general(np.identity(3))
    lat = kwant.lattice.general([(1,0,0),(0,1,0),(0,0,1)])
    sys = kwant.Builder()
    
    sym0 = kwant.TranslationalSymmetry((1, 0, 0))
    lead0 = kwant.Builder(sym0)
    
    sym1 = kwant.TranslationalSymmetry((-1, 0, 0))
    lead1 = kwant.Builder(sym1)
    
    def shape_lead(pos):
        (x, y, z) = pos
        return (0 <= z < H) and (0 <= y < W)
    
    def shape(pos):
        (x, y, z) = pos
        return (0 <= x < L) and (0 <= y < W) and (0 <= z < H)
    
    def hopping_x(site1, site2, par):
        xt, yt, zt = site1.pos
        xs, ys, zs = site2.pos
        return hopx(site1,site2,par) * np.exp(-0.5j * par.Bz * (xt - xs) * (yt + ys))
    
    # scattering system
    sys[lat.shape(shape, (0,0,0))] = lambda site, par: onsite(site, par.scat) - par.mu * np.eye(4)
    sys[kwant.HoppingKind((1,0,0), lat)] = lambda site1, site2, par: hopx(site1, site2, par.scat)
    sys[kwant.HoppingKind((0,1,0), lat)] = lambda site1, site2, par: hopy(site1, site2, par.scat)
    sys[kwant.HoppingKind((0,0,1), lat)] = lambda site1, site2, par: hopz(site1, site2, par.scat)

    # leads
    lead0[lat.shape(shape_lead, (0,0,0))] = lambda site, par: onsite(site, par.lead) - par.mu_lead * np.eye(4)
    lead0[kwant.HoppingKind((1,0,0), lat)] = lambda site1, site2, par: hopx(site1, site2, par.lead)
    lead0[kwant.HoppingKind((0,1,0), lat)] = lambda site1, site2, par: hopy(site1, site2, par.lead)
    lead0[kwant.HoppingKind((0,0,1), lat)] = lambda site1, site2, par: hopz(site1, site2, par.lead)
    
    lead1[lat.shape(shape_lead, (0,0,0))] = lambda site, par: onsite(site, par.lead) - par.mu_lead * np.eye(4)
    lead1[kwant.HoppingKind((1,0,0), lat)] = lambda site1, site2, par: hopx(site1, site2, par.lead)
    lead1[kwant.HoppingKind((0,1,0), lat)] = lambda site1, site2, par: hopy(site1, site2, par.lead)
    lead1[kwant.HoppingKind((0,0,1), lat)] = lambda site1, site2, par: hopz(site1, site2, par.lead)
    
    sys.attach_lead(lead0)
    sys.attach_lead(lead1)
    return sys.finalized() 

Performing the necessary imports.

from __future__ import division, print_function
import numpy as np
import matplotlib
import kwant

import ipywidgets
from IPython.html.widgets import interact
from ipywidgets import StaticInteract, RangeWidget, DropDownWidget
from IPython.display import display_html
from matplotlib import pyplot as plt

import pfaffian as pf
from edx_components import *


:0: FutureWarning: IPython widgets are experimental and may change in the future.

Press this button to show/hide the code used in the notebook:

Table of Contents

In [27]:
Bz = 0.0
mu_lead = 0.7
mu_max = .7

L, W, H = 10, 10, 10
sys = bhz_slab(L, W, H)

momenta = np.linspace(-1.5, 3.5, 60)
par_lead = SimpleNamespace(A1=1.0, A2=1.0, B1=1.0, B2=1.0, C=0.0, D1=0., D2=0., M=1.0, Bz=Bz)
par = SimpleNamespace(lead=par_lead, scat=par_lead, mu=0.0, mu_lead=0.0)
fig, ax = plt.subplots(figsize=((8,6)))
ax.set_color_cycle(['k'])
    
kwant.plotter.bands(sys.leads[0], args=[par], momenta=momenta, show=False, ax=ax);

ax.set_ylim(-3,3);
ax.set_xlim(-2.0,2.0);
ax.set_xlabel('$k$');
ax.set_ylabel('$E$');

vals = np.arange(-1, 2);
ax.set_xticks(vals);
ax.set_xticklabels(["${0}$".format(i) for i in vals]);

#vals = np.arange(-0.8, 1.2, 0.4);
#ax.set_yticks(vals);
#ax.set_yticklabels(["${0}$".format(i) for i in vals]);
In [2]:
Bz = 0.0
mu_lead = 0.7
mu_max = .7

L, W, H = 10, 10, 10
sys = bhz_slab(L, W, H)

momenta = np.linspace(-1.5, 3.5, 60)
par_lead = SimpleNamespace(A1=1.0, A2=1.0, B1=1.0, B2=1.0, C=0.0, D1=0., D2=0., M=-1.0, Bz=Bz)
#par_lead = SimpleNamespace(A1=2.2, A2=4.1, B1=10.0, B2=56.6, C=0.0, D1=1.3, D2=19.6, M=0.28, Bz=Bz)
par = SimpleNamespace(lead=par_lead, scat=par_lead, mu=0.0, mu_lead=0.0)
fig, ax = plt.subplots(figsize=((8,6)))
ax.set_color_cycle(['k'])
    
kwant.plotter.bands(sys.leads[0], args=[par], momenta=momenta, show=False, ax=ax);

ax.set_ylim(-1.5,1.5);
#ax.set_ylim(-40.0,40.0);
#ax.set_xlim(-.6,0.6);
ax.set_xlim(-2.0,2.0);
ax.set_xlabel('$k$');
ax.set_ylabel('$E$');

vals = np.arange(-1, 2);
ax.set_xticks(vals);
ax.set_xticklabels(["${0}$".format(i) for i in vals]);

#vals = np.arange(-0.8, 1.2, 0.4);
#ax.set_yticks(vals);
#ax.set_yticklabels(["${0}$".format(i) for i in vals]);
In [4]:
L, W, H = 10, 10, 10
sys = bhz_slab(L, W, H)
kwant.plot(sys)
Out[4]:
In [5]:
def plot_conductance(sys, energies, args=[par,]):
    # Compute conductance
    data = []
    for energy in energies:
        smatrix = kwant.smatrix(sys, energy, args=[par,])
        data.append(smatrix.transmission(1, 0))

    plt.figure()
    plt.plot(energies, data)
    plt.xlabel("energy [t]")
    plt.ylabel("conductance [e^2/h]")
    plt.show()
In [6]:
Bz = 0.0
mu_lead = 0.7
mu_max = .7
L, W, H = 10, 10, 10
sys = bhz_slab(L, W, H)
par_lead = SimpleNamespace(A1=1.0, A2=1.0, B1=1.0, B2=1.0, C=0.0, D1=0., D2=0., M=-1.0, Bz=Bz)
par = SimpleNamespace(lead=par_lead, scat=par_lead, mu=0.0, mu_lead=0.0)
plot_conductance(sys, energies=[0.01 * i for i in xrange(-50,100)], args=[par,])

/usr/local/lib/python2.7/dist-packages/numpy/core/fromnumeric.py:2499: VisibleDeprecationWarning: `rank` is deprecated; use the `ndim` attribute or function instead. To find the rank of a matrix see `numpy.linalg.matrix_rank`.
  VisibleDeprecationWarning)
In [11]:
Bz = 0.0
mu_lead = 0.7
mu_max = .7
L, W, H = 4, 3, 2
sys = bhz_slab(L, W, H)
par_lead = SimpleNamespace(A1=1.0, A2=1.0, B1=1.0, B2=1.0, C=0.0, D1=0., D2=0., M=-1.0, Bz=Bz)
par = SimpleNamespace(lead=par_lead, scat=par_lead, mu=0.0, mu_lead=0.0)
ham = sys.hamiltonian_submatrix(args=[par,], sparse=True).tocsc()
num_states = 1
energies, states = sl.eigsh(ham, sigma=0, k=num_states)
#print states
#print states.reshape(-1,4, num_states).sum(axis=1)
densities = (np.linalg.norm(states.reshape(-1, 4, num_states), axis=1)**2)
#aaa=kwant.plotter.sys_leads_sites(sys, num_lead_cells=2)
#print aaa
#print densities

fig = plt.figure(figsize=(10, 5))
'''
ax0 = fig.add_subplot(121, projection='3d')
kwant.plot(sys, ax=ax0, site_size=0.5*np.linspace(1, 1, len(densities[:, 0])))
ax0.set_xlim(0, W)
ax0.set_ylim(0, L)
ax0.set_yticks([])
ax0.set_xticks([])
ax0.set_zlim3d([0, H])
ax0.set_zticks([0, H])
ax0.set_zticklabels(['$0$', '$%d$' %H])
'''
densities /= np.max(densities, axis=0, keepdims=True)
print densities[:, 0]
ax1 = fig.add_subplot(projection='3d')
#ax = fig(projection='3d')
kwant.plotter.plot(sys, site_size=densities[:, 0], site_color=densities[:, 0], ax=ax1,
                  cmap='gist_heat_r', colorbar=True, site_lw=0.1)
#kwant.plotter.sys_leads_sites(sys, num_lead_cells=2)
#kwant.plotter.plot(sys, site_size=densities[:, 0], site_color=densities[:, 0], ax=ax,
 #                  cmap='seismic', colorbar=True, site_lw=0, fig_size=(8,5))
#ax1.set_xlim([0, W])
#ax1.set_ylim(0, L)
#ax1.set_yticks([])
#ax1.set_xticks([])
#ax1.set_zlim3d([0, H])
#ax1.set_zticks([0, H])
#ax1.set_zticklabels(['$0$', '$%d$' %H])
##plt.title('Majoranas, $E = $' + scientific_number(abs(energies[0])))
plt.show()

[ 0.26  0.57  1.    1.    0.26  0.46  0.57  0.46  0.46  0.26  0.26  1.
  0.46  0.57  0.26  0.26  0.26  0.26  1.    0.57  0.57  0.57  0.57  0.57]

<matplotlib.figure.Figure at 0x3e80b10>
In [12]:
import tinyarray
from numpy import *
Bz = 0.0
mu_lead = 0.7
mu_max = .7
L, W, H = 4, 4, 3
Nsites = L*W*H
sys = bhz_slab(L, W, H)
par_lead = SimpleNamespace(A1=1.0, A2=1.0, B1=1.0, B2=1.0, C=0.0, D1=0., D2=0., M=-1.0, Bz=Bz)
par = SimpleNamespace(lead=par_lead, scat=par_lead, mu=0.0, mu_lead=0.0)

def density(sys, energy, lead_nr, args=[par,]):
    wf = kwant.wave_function(sys, energy, args=[par,])
    return abs((wf(lead_nr)).reshape(-1,L*W*H).sum(axis=0)**2)

den = density(sys,0.5, 0, args=[par,])
d = zeros((Nsites,))
for i in xrange(Nsites):
    d[i] = den[i]

print d
#kwant.plotter.map(sys, d, cmap='gist_heat_r', colorbar=True)#, site_lw=0)

ldos = kwant.ldos(sys, energy=0.5, args=[par,])
ldos_sum = ldos[::4] + ldos[1::4] + ldos[2::4] + ldos[3::4]
print ldos_sum
#kwant.plotter.plot(sys, local_dos, num_lead_cells=2)

[ 0.82  0.09  0.58  0.18  1.11  0.25  1.    0.11  1.08  0.29  1.01  0.13
  1.19  0.32  0.95  0.33  1.14  0.09  0.71  0.28  0.97  0.07  0.58  0.3
  1.45  0.08  0.91  0.33  0.85  0.03  0.46  0.27  0.69  0.21  0.68  0.04
  0.51  0.06  0.37  0.05  0.8   0.11  0.7   0.07  1.    0.19  0.85  0.13]
[ 0.06  0.05  0.08  0.06  0.03  0.03  0.08  0.03  0.06  0.03  0.03  0.06
  0.06  0.06  0.08  0.05  0.03  0.06  0.08  0.05  0.06  0.03  0.03  0.06
  0.03  0.08  0.03  0.05  0.08  0.06  0.06  0.05  0.03  0.03  0.08  0.05
  0.06  0.06  0.06  0.06  0.08  0.05  0.03  0.03  0.03  0.03  0.06  0.05]
In [ ]: