Hi Hanifa,

The mismatch between the eigenvector size and the number of sites in
the system is because your system has more than one degree of freedom
per site. Please see how to deal with such cases in the tutorial:
https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fkwant-project.org%2Fdoc%2F1%2Ftutorial%2Foperators&data=01%7C01%7Chanifa.tidjani%40kcl.ac.uk%7C4d22e15c6b5845a7382908d729fdcb8f%7C8370cf1416f34c16b83c724071654356%7C0&sdata=bLzi8aLHdVqCfwSycyYE4ZLXhpsJeZR15iUHFF0pZgE%3D&reserved=0

Best,
Anton

On Mon, 26 Aug 2019 at 10:08, Tidjani, Hanifa <hanifa.tidjani@kcl.ac.uk> wrote:
>
>
>
> 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