Wavefunction for 2D NISIN
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
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://kwant-project.org/doc/1/tutorial/operators 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
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 <anton.akhmerov+kd@gmail.com> Sent: Monday, August 26, 2019 9:16:59 AM To: Tidjani, Hanifa <hanifa.tidjani@kcl.ac.uk> Cc: kwant-discuss@kwant-project.org <kwant-discuss@kwant-project.org> 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-proj... 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
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 <hanifa.tidjani@kcl.ac.uk> 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 <anton.akhmerov+kd@gmail.com> Sent: Monday, August 26, 2019 9:16:59 AM To: Tidjani, Hanifa <hanifa.tidjani@kcl.ac.uk> Cc: kwant-discuss@kwant-project.org <kwant-discuss@kwant-project.org> 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-proj...
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
Hi Anton, I still don’t understand exactly how I can implement this. From what I can understand you choose the energy not the eigenstate for the density function. Is there a way to use the eigenvector solution in 2.4 for the wavefunction for a system with two sublattices? Or am I simply misunderstanding where in the density function you encode the eigenstate? Kind regards, Hanifa From: Anton Akhmerov <anton.akhmerov+kd@gmail.com> Sent: Thursday, August 29, 2019 2:04:47 PM To: Tidjani, Hanifa <hanifa.tidjani@kcl.ac.uk> Cc: kwant-discuss@kwant-project.org <kwant-discuss@kwant-project.org> Subject: Re: [Kwant] Wavefunction for 2D NISIN 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 <hanifa.tidjani@kcl.ac.uk> 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 <anton.akhmerov+kd@gmail.com> Sent: Monday, August 26, 2019 9:16:59 AM To: Tidjani, Hanifa <hanifa.tidjani@kcl.ac.uk> Cc: kwant-discuss@kwant-project.org <kwant-discuss@kwant-project.org> 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%7Cfe0e226f1cb240cc912608d72c817fe2%7C8370cf1416f34c16b83c724071654356%7C0&sdata=Nr9oHbH%2FRzVmWJaYBwpYapAaXO3I%2BZZTIR6ui3yWw0g%3D&reserved=0<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%7Cfe0e226f1cb240cc912608d72c817fe2%7C8370cf1416f34c16b83c724071654356%7C0&sdata=Nr9oHbH%2FRzVmWJaYBwpYapAaXO3I%2BZZTIR6ui3yWw0g%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
Kind regards, Hanifa.
Hi,
I still don’t understand exactly how I can implement this. From what I can understand you choose the energy not the eigenstate for the density function.
I am unclear what you mean here. In your code you do the following: ham_mat = sys.hamiltonian_submatrix(sparse=True, params = params) evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0)) kwant.plotter.map(sys, np.abs(evecs[:, 9])**2 You take the hamiltonian of the scattering region, diagonalize it, and then plot *a single eigenvector*. Your problem was that because you have
1 degree of freedom per site you cannot simply plot the eigenvector as-is. You need to first compute the density per site for this eigenvector, which you can do using the operator module:
rho = kwant.operator.Density(np.kron(-s_z, s_0)) # holes count +1 to the density, and electrons -1 kwant.plotter.map(sys, rho(evecs[:, 9)) Note that you call the density operator with a single eigenvector. Independently of the above it's not clear to me that you should be diagonalizing just the scattering region Hamiltonian, given that your system has leads. It seems to me that you should be calculating the scattering states and the density of those. Happy Kwanting, Joe
participants (3)
-
Anton Akhmerov -
Joseph Weston -
Tidjani, Hanifa