Mailman 3 October 2020 - Kwant-discuss

Implementing the BHZ model
by Victor Regis Oct. 12, 2020

Oct. 12, 2020

Hello all, I'm trying to implement the BHZ model already in the tight-biding model according to the paper [Phys. Rev. B 96, 155302 (2017) or arXiv:1707.03624]. I use the equations (5),(6) and (7), however with no magnetic field. I'm using Kwant's documentation example 2.10 as a guideline (For what I have to do next it's more convenient that I don't use the discretize method). After some struggling I got the band structure of both mine and Kwant's example to be the same, however my spin density is reversed, it's up where it's supposed to be down. So far I haven't figured out why. Also, the main difference between the implementation is only the make_syst part, all the plotting in my code are snippets from Kwant's discretize.py; My code is below: import kwant import matplotlib.pyplot as plt import numpy as np import random as rd from math import * import matplotlib as mpl import tinyarray sigma_0 = tinyarray.array([[1, 0], [0, 1]]) sigma_x = tinyarray.array([[0, 1], [1, 0]]) sigma_y = tinyarray.array([[0, -1j], [1j, 0]]) sigma_z = tinyarray.array([[1, 0], [0, -1]]) Id = np.eye(4) G1 = np.kron(sigma_z,sigma_x) G2 = -np.kron(sigma_0,sigma_y) G5 = np.kron(sigma_0,sigma_z) def make_syst(a=20,A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0,L=200,W=200): lat = kwant.lattice.square(a,norbs=4) syst = kwant.Builder() def square(pos): (x, y) = pos return (0 <= y < W and 0 <= x < L) def onsite(site): (x, y) = site.pos return C*Id+M*G5-(4/a**2)*(D*Id+B*G5) #equation 5 def hoppingx(site1,site2): return (1/a**2)*(D*Id+B*G5)+(1j/(2*a))*(A*G1) #eq 6 def hoppingy(site1,site2): return (1/a**2)*(D*Id+B*G5)-(1j/(2*a))*(A*G2) #eq 7 syst[lat.shape(square,(0,0))] = onsite syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hoppingx syst[kwant.builder.HoppingKind((0,1), lat, lat)] = hoppingy def lead_shape(pos): (x, y) = pos return (0 <= y < W) #### Define the left lead. #### lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0))) lead[lat.shape(lead_shape,(0,0))] = onsite # hoppings in x-direction lead[kwant.builder.HoppingKind((1, 0), lat, lat)] =hoppingx # hoppings in y-directions lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = hoppingy return syst,[lead,lead.reversed()] def main(): syst,leads = make_syst(W=1000,L=2000) # Attach the leads to the system. for lead in leads: syst.attach_lead(lead) # Finalize the system. syst = syst.finalized() # Check that the system looks as intended. #kwant.plot(syst) # Compute the band structure of lead 0. kwant.plotter.bands(syst.leads[0], momenta=np.linspace(-0.3, 0.3, 201), show=False) plt.grid() plt.xlim(-.3, 0.3) plt.ylim(-0.05, 0.05) plt.xlabel('momentum [1/A]') plt.ylabel('energy [eV]') plt.show() wf = kwant.wave_function(syst, energy=0) # prepare density operators sigma_z = np.array([[1, 0], [0, -1]]) prob_density = kwant.operator.Density(syst, np.eye(4)) spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2))) # calculate expectation values and plot them wf_sqr = sum(prob_density(psi) for psi in wf(0)) # states from left lead rho_sz = sum(spin_density(psi) for psi in wf(0)) # states from left lead fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(16, 4)) kwant.plotter.map(syst, wf_sqr, ax=ax1) kwant.plotter.map(syst, rho_sz, ax=ax2) ax = ax1 im = [obj for obj in ax.get_children() if isinstance(obj, mpl.image.AxesImage)][0] fig.colorbar(im, ax=ax) ax = ax2 im = [obj for obj in ax.get_children() if isinstance(obj, mpl.image.AxesImage)][0] fig.colorbar(im, ax=ax) ax1.set_title('Probability density') ax2.set_title('Spin density') plt.show() # Call the main function if the script gets executed (as opposed to imported). # See <http://docs.python.org/library/__main__.html>. if __name__ == '__main__': main() Regards, Victor Regis

1 0

Implementing the BHZ model
by Victor Regis Oct. 12, 2020

Oct. 12, 2020

Hello all, I'm trying to implement the BHZ model already in the tight-biding model according to the paper [Phys. Rev. B 96, 155302 (2017) or arXiv:1707.03624]. I use the equations (5),(6) and (7), however with no magnetic field. I'm using Kwant's documentation example 2.10 as a guideline (For what I have to do next it's more convenient that I don't use the discretize method). After some struggling I got the band structure of both mine and Kwant's example to be the same, however my spin density is reversed, it's up where it's supposed to be down. So far I haven't figured out why. Also, the main difference between the implementation is only the make_syst part, all the plotting in my code are snippets from Kwant's discretize.py; My code is below: import kwant import matplotlib.pyplot as plt import numpy as np import random as rd from math import * import matplotlib as mpl import tinyarray sigma_0 = tinyarray.array([[1, 0], [0, 1]]) sigma_x = tinyarray.array([[0, 1], [1, 0]]) sigma_y = tinyarray.array([[0, -1j], [1j, 0]]) sigma_z = tinyarray.array([[1, 0], [0, -1]]) Id = np.eye(4) G1 = np.kron(sigma_z,sigma_x) G2 = -np.kron(sigma_0,sigma_y) G5 = np.kron(sigma_0,sigma_z) def make_syst(a=20,A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0,L=200,W=200): lat = kwant.lattice.square(a,norbs=4) syst = kwant.Builder() def square(pos): (x, y) = pos return (0 <= y < W and 0 <= x < L) def onsite(site): (x, y) = site.pos return C*Id+M*G5-(4/a**2)*(D*Id+B*G5) #equation 5 def hoppingx(site1,site2): return (1/a**2)*(D*Id+B*G5)+(1j/(2*a))*(A*G1) #eq 6 def hoppingy(site1,site2): return (1/a**2)*(D*Id+B*G5)-(1j/(2*a))*(A*G2) #eq 7 syst[lat.shape(square,(0,0))] = onsite syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hoppingx syst[kwant.builder.HoppingKind((0,1), lat, lat)] = hoppingy def lead_shape(pos): (x, y) = pos return (0 <= y < W) #### Define the left lead. #### lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0))) lead[lat.shape(lead_shape,(0,0))] = onsite # hoppings in x-direction lead[kwant.builder.HoppingKind((1, 0), lat, lat)] =hoppingx # hoppings in y-directions lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = hoppingy return syst,[lead,lead.reversed()] def main(): syst,leads = make_syst(W=1000,L=2000) # Attach the leads to the system. for lead in leads: syst.attach_lead(lead) # Finalize the system. syst = syst.finalized() # Check that the system looks as intended. #kwant.plot(syst) # Compute the band structure of lead 0. kwant.plotter.bands(syst.leads[0], momenta=np.linspace(-0.3, 0.3, 201), show=False) plt.grid() plt.xlim(-.3, 0.3) plt.ylim(-0.05, 0.05) plt.xlabel('momentum [1/A]') plt.ylabel('energy [eV]') plt.show() wf = kwant.wave_function(syst, energy=0) # prepare density operators sigma_z = np.array([[1, 0], [0, -1]]) prob_density = kwant.operator.Density(syst, np.eye(4)) spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2))) # calculate expectation values and plot them wf_sqr = sum(prob_density(psi) for psi in wf(0)) # states from left lead rho_sz = sum(spin_density(psi) for psi in wf(0)) # states from left lead fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(16, 4)) kwant.plotter.map(syst, wf_sqr, ax=ax1) kwant.plotter.map(syst, rho_sz, ax=ax2) ax = ax1 im = [obj for obj in ax.get_children() if isinstance(obj, mpl.image.AxesImage)][0] fig.colorbar(im, ax=ax) ax = ax2 im = [obj for obj in ax.get_children() if isinstance(obj, mpl.image.AxesImage)][0] fig.colorbar(im, ax=ax) ax1.set_title('Probability density') ax2.set_title('Spin density') plt.show() # Call the main function if the script gets executed (as opposed to imported). # See <http://docs.python.org/library/__main__.html>. if __name__ == '__main__': main() Regards, Victor Regis

1 0

Plotting bands using wraparound for 3D system
by Naveen Yadav Oct. 2, 2020

Oct. 2, 2020

Dear KWANT developers, I am trying to plot the energy bands as a function of k_x. The system has finite width in z direction and infinite in y and x-direction. I have used the wraparound module to plot the bands and I kept the translation invariance in x-direction and attached the leads. But while plotting the bands of lead using kwant.plotter.bands it shows an error. Please suggest to me the correct way to plot the bands for such a system. Thank you. Regards import numpy from kwant.wraparound import wraparound,plot_2d_bands from matplotlib import pyplot import pandas as pd get_ipython().magic('matplotlib inline') import sys import itertools import collections import cmath import numpy as np import tinyarray from numpy import pi,sqrt,arccos,array, exp import kwant from kwant.builder import herm_conj import scipy.linalg import scipy.sparse.linalg as sla import matplotlib as mpl mpl.rc('text', usetex = True) mpl.rc('font', family = 'serif') #%matplotlib notebook sigma_0 = tinyarray.array([[1, 0], [0, 1]]) sigma_x = tinyarray.array([[0, 1], [1, 0]]) sigma_y = tinyarray.array([[0, -1j], [1j, 0]]) sigma_z = tinyarray.array([[1, 0], [0, -1]]) a=1; L=30; W=1; H=1; def onsite(site): return (2.0 * np.cos(1.5) + 4) * sigma_z def hoppingx(site0, site1): return (-sigma_z - 1j * sigma_x) def hoppingy(site0, site1): return -sigma_z - 1j * sigma_y def hoppingz(site0, site1): return (-1 * 1.0 * sigma_z) lat =kwant.lattice.cubic(norbs=2) syst = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((0, W, 0)))) syst[(lat(x, y, z) for x in range(L) for y in range(W) for z in range(H))] = onsite syst[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingx syst[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy syst[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingz syst = kwant.wraparound.wraparound(syst) #Left lead leadL = kwant.Builder(kwant.TranslationalSymmetry((-a, 0, 0), lat.vec((0, W, 0)))) #trans inv in transport dir and ky leadL[(lat(x, y, z) for x in range(L) for y in range(W) for z in range(H))] = onsite leadL[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingx leadL[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy leadL[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingz #Right lead leadR = kwant.Builder(kwant.TranslationalSymmetry((-a, 0, 0), lat.vec((0, W, 0)))) leadR[(lat(x, y, z) for x in range(L) for y in range(W) for z in range(H))] = onsite leadR[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingx leadR[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy leadR[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingz leadL = kwant.wraparound.wraparound(leadL, keep=0) #keep the translational invariance in x-direction leadR = kwant.wraparound.wraparound(leadR, keep=0) syst.attach_lead(leadL) syst.attach_lead(leadR) kwant.plot(syst) syst = syst.finalized() kwant.plotter.bands(syst.leads[1], momenta=np.linspace(-pi, pi, 101, endpoint=False)) -

1 0