I am trying to plot the wavefunction for a 2D NISIN Hamiltonian however I am unable to get a result, this is my code:
import kwant
import tinyarray as ta
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
import winsound
import os
from datetime import datetime
import scipy.sparse.linalg as sla
sig_0 = np.identity(4)
s_0 = np.identity(2)
s_z = np.array([[1., 0.], [0., -1.]])
s_x = np.array([[0., 1.], [1., 0.]])
s_y = np.array([[0., -1j], [1j, 0.]])
tau_z = ta.array(np.kron(s_z, s_0))
tau_x = ta.array(np.kron(s_x, s_0))
tau_y = ta.array(np.kron(s_y, s_0))
sig_z = ta.array(np.kron(s_0, s_z))
tau_zsig_x = ta.array(np.kron(s_z, s_x))
tau_zsig_y = ta.array(np.kron(s_z, s_y))
# Lorentzian definition
def lorentzianxy(x, y, ind,y0):
gam = 3
L = params["L"]
#if 0<ind<L:
return gam ** 2 /( gam ** 2 + (x - ind) ** 2 + (y - y0) ** 2)
# else :
# return 0
# Define potential
def potential(site, pot,ind,y0):
(x, y) = site.pos
return pot * lorentzianxy(x,y,ind,y0)
# Make system
def makeNISIN2D(params):
# List all the params
W=params["W"]
L=params["L"]
pot = params["pot"]
ind = params["ind"]
mu=params["mu"]
B=params["B"]
Delta=params["Delta"]
alpha=params["alpha"]
t=params["t"]
barrier = params["barrier"]
Smu = params["Smu"]
def lorentzianxy(x, y, ind,y0):
gam = 3
L = params["L"]
#if 0<ind<L:
return gam ** 2 /( gam ** 2 + (x - ind) ** 2 + (y - y0) ** 2)
# else :
# return 0
# Define potential
def potential(site, pot,ind,y0):
(x, y) = site.pos
return pot * lorentzianxy(x,y,ind,y0)
# if y < lorentzian(x,ind)*4:
#
# return (tru_pot)
# else:
# return 0
#
def Straightpotential(site, pot,ind):
(x, y) = site.pos
if x == ind:
return pot
else:
return 0
def onsiteSc(site, pot,ind,y0):
#(x,y)=site.pos
#L=params["L"]
# if x>L/2:
return (4 * t - Smu) * tau_z + B * sig_z + Delta * tau_x - potential(site, pot,ind,y0)*tau_z
# else:
#return (4 * t ) * tau_z + B * sig_z + Delta * tau_x - potential(site, pot,ind,y0)*tau_z
def onsiteNormal(site, pot,ind,y0):
return (4 * t - mu) * tau_z #- potential(site, pot,ind,y0)*tau_z
def onsiteBarrier(site, pot,ind,y0):
return (4 * t - mu + barrier) * tau_z #- potential(site, pot,ind,y0)*tau_z
def hop_x(site, pot, ind,y0):
return -t * tau_z + 0.5j * alpha * tau_zsig_y
def hop_y(site, pot, ind,y0):
return -t * tau_z - .5j * alpha * tau_zsig_x
# On each site, electron and hole orbitals.
lat = kwant.lattice.square(norbs=4)
syst = kwant.Builder()
# S
syst[(lat(i, j) for i in range(1,L-1) for j in range(W))] = onsiteSc
barrierLen = 1;
# I's
syst[(lat(i, j) for i in range(barrierLen) for j in range(W))] = onsiteBarrier
syst[(lat(i, j) for i in range(L-barrierLen, L)for j in range(W))] = onsiteBarrier
syst[kwant.builder.HoppingKind((1, 0), lat, lat)]= hop_x
syst[kwant.builder.HoppingKind((0, 1), lat, lat)]= hop_y
# N's
lead = kwant.Builder(kwant.TranslationalSymmetry((-1,0)),
conservation_law=-tau_z#, particle_hole=tau_y
)
lead[(lat(0, j) for j in range(W))] = onsiteNormal
lead[kwant.builder.HoppingKind((1, 0), lat, lat)]= hop_x
lead[kwant.builder.HoppingKind((0, 1), lat, lat)]= hop_y
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
return syst
def sorted_eigs(ev):
evals, evecs = ev
evals, evecs = map(np.array, zip(*sorted(zip(evals, evecs.transpose()))))
return evals, evecs.transpose()
params = dict(mu=0.3, Delta=.1, alpha=.8, t=1.0, barrier = 2.0, pot = 0.0,W = 5, L = 30, ind = 0, B = 0.3, Smu=0.00,y0=0)
sys = makeNISIN2D(params)
sys = sys.finalized()
# Calculate the wave functions in the system.
ham_mat = sys.hamiltonian_submatrix(sparse=True, params = params)
evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))
# Plot the probability density of the 10th eigenmode.
kwant.plotter.map(sys, np.abs(evecs[:, 9])**2,
colorbar=False, oversampling=1)
And the error I get is this: "ValueError: The number of sites doesn't match the number of provided values."
When I print the values they indeed are more than what they should be. I've tried a couple of things but I'm quite at a loss of what to do. I want to plot the wavefunction so I can see the Majoranas at the edges of the superconductor so that I can see how they respond to a local perturbation in the potential. Any tips would be very helpful.
Kind regards,
Hanifa
