Hi Hanifa, This is certainly possible. It's up to you for which wave function you compute the density. You can compute the eigenstates and compute the expectation all the same.
Best, Anton On Mon, 26 Aug 2019 at 18:46, Tidjani, Hanifa <[email protected]> wrote: > > Dear Anton, > > > > Thank you for your reply. > > > > This does indeed work, however I am concerned that I am driving it at a > certain energy and then plotting the wave function. Is there a way to plot > the wave function of a particular eigenstate rather than at a particular > energy as is done in tutorial 2.4? > > > > Kind regards, > > > > Hanifa > > > > ________________________________ > From: Anton Akhmerov <[email protected]> > Sent: Monday, August 26, 2019 9:16:59 AM > To: Tidjani, Hanifa <[email protected]> > Cc: [email protected] <[email protected]> > Subject: Re: [Kwant] Wavefunction for 2D NISIN > > 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 <[email protected]> > 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
