Mailman 3 June 2023 - Kwant-discuss

Package announcement: Pymablock
by Anton Akhmerov 12 Nov '23

12 Nov '23

Hello fellow Kwanters, I am happy to announce Pymablock, a Python package useful to combine with Kwant. Its main purpose is constructing effective models using perturbation theory, which may help you in different ways. Pymablock allows you to analyse the behavior of a few levels of a large Hamiltonian, or to construct a symbolic k.p model. We originally wrote the code for one of our projects, but then saw that this tool is useful almost for the entire group, so we decided to turn it into a package. We spent countless hours generalizing, optimizing, and testing it so that it is efficient and reliable. Now it handles multiple perturbations to any order, and it is compatible with symbolic and numerical calculations. We hope that Pymablock will save you a lot of effort and computation time. To see where it can help, check out our its documentation (https://pymablock.rtfd.io/), and in particular the tutorials about the k.p model of bilayer graphene (https://pymablock.readthedocs.io/en/latest/tutorial/bilayer_graphene.html) or the induced superconducting gap (https://pymablock.readthedocs.io/en/latest/tutorial/induced_gap.html). Just like Kwant, Pymablock is open source (https://gitlab.kwant-project.org/qt/pymablock), and you are welcome to follow and contribute to its development. Let us know if you have any questions or comments! Also if there's interest we can organize an online introduction of the package Best, Anton Akhmerov and the Pymablock maintainers: Isidora Araya Day and Sebastian Miles

3 3

Realistic values of hopping parameters
by Yu Li 18 Oct '23

18 Oct '23

Dear all, Now I’m studying some transport property based on multilayered system, with nonmagnetic layer/ferromagnetic layer (NM/FM) stack. May I ask when creating a tight binding system, and setting hoppings by the following function: syst[(lat(1, 0), lat(1, 1))] = t, how to associate the hopping with a realistic system? Such as hoppings between FM and FM atoms, or between FM and NM atoms? I found in many literatures people just simply put t=1, but I assume the value should depend on the type of atoms. My another question is, is it possible to define more than one type of hopping between two sites? For the tight-binding approximation of p electrons, there are two types of hopping integrals: Vppσ, and Vppπ, so I’m wondering if it’s possible to simulate the effect of both hoppings? Thanks a lot for reading my questions! Cheers, Yu Li

4 3

Scattering region and lead with different kind of lattice
by Gabriel Garcia 11 Jul '23

11 Jul '23

Hi devs, I am learning kwant alone, but i have doubts. I am trying to construct a system with a diatomic chain as scattering region and the leads as a monoatomic chain. But i am doing something wrong. Someone else can help me? Bellow i will show the code: ##################################################################################################################### from matplotlib import pyplot as plt from matplotlib import backends import kwant # Defining the parameters L = 5 # range of scattering region C = 2 # primitive vector # Definindo o parametro de Hopping t = 5 # hopping parameter lat = kwant.lattice.general([(C,0),(0,1)],[(0,0),(1,0)], norbs=1) a, b = lat.sublattices syst = kwant.Builder() # Scattering region for i in range(L): syst[a(i, 0)] = 0 syst[b(i, 0)] = 0 syst[kwant.builder.HoppingKind((0, 0), a, b)] = -t syst[kwant.builder.HoppingKind((1, 0), b, a)] = -t kwant.plot(syst) plt.show() # leaft lead (lead 1) lat_lead1 = kwant.lattice.chain(C, norbs=1) sym_lead1 = kwant.TranslationalSymmetry((-C, 0)) lead1 = kwant.Builder(sym_lead1) lead1[lat_lead(0,0)] = 4*t lead1[lat_lead1.neighbors()] = -t syst[(lat_lead1(-C, 0), a(0, 0))] = -t syst.attach_lead(lead1) kwant.plot(syst) --------------------------------------------------------------------------- ValueError Traceback (most recent call last) ~\anaconda3\lib\site-packages\kwant\builder.py in __new__(cls, family, tag, _i_know_what_i_do) 75 try: ---> 76 tag = family.normalize_tag(tag) 77 except (TypeError, ValueError) as e: ~\anaconda3\lib\site-packages\kwant\lattice.py in normalize_tag(self, tag) 475 if len(tag) != self.lattice_dim: --> 476 raise ValueError("Dimensionality mismatch.") 477 return tag ValueError: Dimensionality mismatch. During handling of the above exception, another exception occurred: ValueError Traceback (most recent call last) <ipython-input-46-3d36b05d2870> in <module> 8 lead1[lat_lead1.neighbors()] = -t 9 ---> 10 syst[(lat_lead1(-C, 0), a(0, 0))] = -t 11 12 syst.attach_lead(lead1) ~\anaconda3\lib\site-packages\kwant\builder.py in __call__(self, *tag) 195 raise ValueError('Use site_family(1, 2) instead of ' 196 'site_family((1, 2))!') --> 197 return Site(self, tag) 198 199 ~\anaconda3\lib\site-packages\kwant\builder.py in __new__(cls, family, tag, _i_know_what_i_do) 77 except (TypeError, ValueError) as e: 78 msg = 'Tag {0} is not allowed for site family {1}: {2}' ---> 79 raise type(e)(msg.format(repr(tag), repr(family), e.args[0])) 80 return tuple.__new__(cls, (family, tag)) 81 ValueError: Tag (-2, 0) is not allowed for site family kwant.lattice.Monatomic([[2.0]], [0.0], '', 1): Dimensionality mismatch. ############################################################################################# Thanks! Gabriel Garcia

2 8

Plotting Particle-Hole density in KWANT
by sc20rs001＠iiserkol.ac.in 30 Jun '23

30 Jun '23

Hi KWANT users! I was trying to plot particle-hole density in KWANT. My system has spin, pseudospin and particle-hole so, norbs = 8. I used kwant.operator.Density(np.kron(np.eye(4), \tau_z)) to plot particle-hole density but could not obtain correct result. It would be great if someone could tell how to plot particle hole density in this case ? Thanks in advance.

1 0

random sites
by Hossein Karbaschi 20 Jun '23

20 Jun '23

Hi everybody, In the following code I want to delete few random sites from the scattering region. Please help me in this case. Best regards, Hossein import kwant from math import sqrt import matplotlib.pyplot as plt import tinyarray import numpy as np import math import matplotlib #Constants......................................................... #Parameter......................................................... a0=2.617*sqrt(3); d0=2.617; dz0=1.59047; D0=sqrt((d0)**2+(dz0)**2); #st=0; st=-0.0; dz=1.59047+st; d=-0.3*dz+3.09; a1=-0.3*sqrt(3)*dz+5.35; D=sqrt((d)**2+(dz)**2); scale=(D/D0)**(-8); qo=[]; N=dz/D; delta1=0j; delta2=0; R=(2*a1)**2; #on-site energy................................................... Es=-10.906; Ep=-0.486; Ep1=-0.486; Ep2=-0.486; #nearest neighbor................................................ sssigma=-0.608*scale; spsigma=1.32*scale; ppsigma=1.854*scale; pppay=-0.6*scale; #next nearest neighbor......................................... ppsigma2=0.156*scale; pppay2=0; sssigma2=0; spsigma2=0; #onsite potentials ........................................... pot_a=10; pot_b=0; #nearest neighbors............................................................. l1=(math.cos(math.pi/6))*(d/D); m1=(math.sin(math.pi/6))*(d/D); N1=dz/D; #second neighbors.............................................................. l2=(math.cos(math.pi/3)) m2=(math.sin(math.pi/3)) latt = kwant.lattice.general([(a1,0),(a1*0.5,a1*math.sqrt(3)/2)], [(a1/2,-d/2),(a1/2,d/2)]) a,b = latt.sublattices syst= kwant.Builder() # Onsite matrice energt Onsite_a=tinyarray.array([[Es,0,0,0],[0,Ep,delta1,0],[0,-delta1,Ep,0], [0,0,0,Ep1]]) Onsite_b=tinyarray.array([[Es,0,0,0],[0,Ep,delta1,0],[0,-delta1,Ep,0], [0,0,0,Ep2]]) pot_edge_a=tinyarray.array([[Es+pot_a,0,0,0],[0,Ep+pot_a,delta1,0],[0,-delta1,Ep+pot_a,0], [0,0,0,Ep1+pot_a]]) pot_edge_b=tinyarray.array([[Es+pot_b,0,0,0],[0,Ep+pot_b,delta1,0],[0,-delta1,Ep+pot_b,0], [0,0,0,Ep2+pot_b]]) #nearest neighbors............................................................. first_up=tinyarray.array([[sssigma, l1*spsigma,m1*spsigma,N1*spsigma], [-l1*spsigma, l1**2*ppsigma+(1-l1**2)*pppay,l1*m1*(ppsigma-pppay),l1*N1*(ppsigma-pppay)], [-m1*spsigma ,l1*m1*(ppsigma-pppay), m1**2*ppsigma+(1-m1**2)*pppay,N1*m1*(ppsigma-pppay)], [-N1*spsigma,l1*N1*(ppsigma-pppay),N1*m1*(ppsigma-pppay),N1**2*ppsigma+(1-N1**2)*pppay]]) first_left_up=tinyarray.array([[sssigma,-l1*spsigma,m1*spsigma,N1*spsigma], [l1*spsigma, l1**2*ppsigma+(1-l1**2)*pppay,- l1*m1*(ppsigma-pppay),-l1*N1*(ppsigma-pppay)], [-m1*spsigma ,-l1*m1*(ppsigma-pppay), m1**2*ppsigma+(1-m1**2)*pppay,N1*m1*(ppsigma-pppay)], [-N1*spsigma,-l1*N1*(ppsigma-pppay),N1*m1*(ppsigma-pppay),N1**2*ppsigma+(1-N1**2)*pppay]]) first=tinyarray.array([[sssigma,0,-(d/D)*spsigma,N1*spsigma], [0,pppay,0,0], [(d/D)*spsigma ,0, (d/D)**2*ppsigma+(1-(d/D)**2)*pppay,-N1*(d/D)*(ppsigma-pppay)], [-N1*spsigma,0,-N1*(d/D)*(ppsigma-pppay),N1**2*ppsigma+(1-N1**2)*pppay]]) #second neighbors.............................................................. second_up=tinyarray.array([[sssigma2,l2*spsigma2,m2*spsigma2,0], [-l2*spsigma2,l2**2*ppsigma2+(1-l2**2)*pppay2,l2*m2*(ppsigma2-pppay2),0], [-m2*spsigma2 ,l2*m2*(ppsigma2-pppay2),m2**2*ppsigma2+(1-m2**2)*pppay2,0], [0,0,0,pppay2]]) second_left_up=tinyarray.array([[sssigma2,-l2*spsigma2,m2*spsigma2,0], [l2*spsigma2,l2**2*ppsigma2+(1-l2**2)*pppay2,-l2*m2*(ppsigma2-pppay2),0], [-m2*spsigma2 ,-l2*m2*(ppsigma2-pppay2),m2**2*ppsigma2+(1-m2**2)*pppay2,0], [0,0,0,pppay2]]) second_right_down=tinyarray.array([[sssigma2,l2*spsigma2,-m2*spsigma2,0], [-l2*spsigma2,l2**2*ppsigma2+(1-l2**2)*pppay2,-l2*m2*(ppsigma2-pppay2),0], [m2*spsigma2,-l2*m2*(ppsigma2-pppay2),m2**2*ppsigma2+(1-m2**2)*pppay2,0], [0,0,0,pppay2]]) second_down=tinyarray.array([[sssigma2,-l2*spsigma2,-m2*spsigma2,0], [l2*spsigma2,l2**2*ppsigma2+(1-l2**2)*pppay2,l2*m2*(ppsigma2-pppay2),0], [m2*spsigma2 ,l2*m2*(ppsigma2-pppay2),m2**2*ppsigma2+(1-m2**2)*pppay2,0], [0,0,0,pppay2]]) second_right=tinyarray.array([[sssigma2,spsigma2,0,0], [-spsigma2,ppsigma2,0,0],[0 ,0,pppay2,0], [0,0,0,pppay2]]) second_left=tinyarray.array([[sssigma2,-spsigma2,0,0], [spsigma2,ppsigma2,0,0],[0 ,0,pppay2,0], [0,0,0,pppay2]]) #................................................................................... def rectangle(pos): x, y = pos # z=x**2+y**2 return -9.2*a1<x<9.2*a1 and -12*d<=y<12*d syst[a.shape(rectangle, (1,1))] = Onsite_a syst[b.shape(rectangle, (1,1))] = Onsite_b #nearest neighbors............................................................. syst[[kwant.builder.HoppingKind((0,0),a,b)]] = first syst[[kwant.builder.HoppingKind((0,1),a,b)]] = first_up syst[[kwant.builder.HoppingKind((-1,1),a,b)]] = first_left_up #second neighbors.............................................................. syst[[kwant.builder.HoppingKind((0,1),a,a)]]=second_up syst[[kwant.builder.HoppingKind((-1,1),a,a)]]=second_left_up syst[[kwant.builder.HoppingKind((1,-1),a,a)]]=second_right_down syst[[kwant.builder.HoppingKind((0,-1),a,a)]]=second_down syst[[kwant.builder.HoppingKind((1,0),a,a)]]=second_right syst[[kwant.builder.HoppingKind((-1,0),a,a)]]=second_left syst[[kwant.builder.HoppingKind((0,1),b,b)]]=second_up syst[[kwant.builder.HoppingKind((-1,1),b,b)]]=second_left_up syst[[kwant.builder.HoppingKind((1,-1),b,b)]]=second_right_down syst[[kwant.builder.HoppingKind((0,-1),b,b)]]=second_down syst[[kwant.builder.HoppingKind((1,0),b,b)]]=second_right syst[[kwant.builder.HoppingKind((-1,0),b,b)]]=second_left sym = kwant.TranslationalSymmetry(latt.vec((1,0))) sym.add_site_family(latt.sublattices[0], other_vectors=[(-1, 2)]) sym.add_site_family(latt.sublattices[1], other_vectors=[(-1, 2)]) lead = kwant.Builder(sym) # ---- random impurity # --- leads def lead_shape(pos): x, y = pos return -12*d<=y<12*d and -5.2*a1<x<5.2*a1 lead[a.shape(lead_shape, (1,1))] = Onsite_a lead[b.shape(lead_shape, (1,1))] = Onsite_b lead[[kwant.builder.HoppingKind((0,0),a,b)]] =first lead[[kwant.builder.HoppingKind((0,1),a,b)]] =first_up lead[[kwant.builder.HoppingKind((-1,1),a,b)]] =first_left_up lead[[kwant.builder.HoppingKind((0,1),a,a)]]=second_up lead[[kwant.builder.HoppingKind((-1,1),a,a)]]=second_left_up lead[[kwant.builder.HoppingKind((1,-1),a,a)]]=second_right_down lead[[kwant.builder.HoppingKind((0,-1),a,a)]]=second_down lead[[kwant.builder.HoppingKind((1,0),a,a)]]=second_right lead[[kwant.builder.HoppingKind((-1,0),a,a)]]=second_left lead[[kwant.builder.HoppingKind((0,1),b,b)]]=second_up lead[[kwant.builder.HoppingKind((-1,1),b,b)]]=second_left_up lead[[kwant.builder.HoppingKind((1,-1),b,b)]]=second_right_down lead[[kwant.builder.HoppingKind((0,-1),b,b)]]=second_down lead[[kwant.builder.HoppingKind((1,0),b,b)]]=second_right lead[[kwant.builder.HoppingKind((-1,0),b,b)]]=second_left syst.attach_lead(lead,add_cells=0) syst.attach_lead(lead.reversed(),add_cells=0) fsys = syst.finalized() def family_color(site): #print(sys[site]) if syst[site]==Onsite_a: return 'black' elif syst[site]==Onsite_b: return 'green' elif syst[site]==pot_edge_a: return 'yellow' else: return 'blue' ax=kwant.plot(syst,site_color=family_color); ax.savefig('syst.pdf')

2 1

Regarding Conductance in Graphene ( Andreev Reflection)
by Sarbajit Mazumdar 15 Jun '23

15 Jun '23

Hi All, I am trying to reproduce Fig. 4 in PRL 97, 067007 (2006) - Specular Andreev Reflection in Graphene. This had been done in KWANT in the PhD thesis of Tibor Sekera - "Quantum transport of fermions in honeycomb lattices and cold atomic systems", chapter 4. My code is the following. Can Anyone help me? I am not getting same result. I am trying the same code this question has been asked before here (https://mail.python.org/archives/list/kwant-discuss@python.org/thread/3JH5W…), but no fruitful outcome was gained. I shall be highly grateful if anyone can help me. ################################################################ import kwant, tinyarray from math import pi, sqrt from matplotlib import pyplot as plt import numpy as np plt.rcParams['figure.dpi']=150 ####___Fundamental Constants___#### e = 1.60217662 * (10**-19) hbar = 1.0545718 * (10**-34) h = 6.62607004 * (10**-34) e_hb = e / hbar ################################### # Specifying 2x2 Pauli matrices acting on spin s_x = tinyarray.array([[0, 1], [1, 0]]) s_y = tinyarray.array([[0, -1j], [1j, 0]]) s_z = tinyarray.array([[1, 0], [0, -1]]) s_0 = tinyarray.array([[1, 0], [0, 1]]) # Specifying 2x2 Pauli matrices acting on particle-hole tau_0 = tinyarray.array([[1, 0], [0, 1]]) tau_x = tinyarray.array([[0, 1], [1, 0]]) tau_y = tinyarray.array([[0, -1j], [1j, 0]]) tau_z = tinyarray.array([[1, 0], [0, -1]]) #Fix hopping, Fermi energy, and superconducting gap t = 1 # Ef = 0 Ef_S=2 Delta= 2 a=1 V0=0.5 #Length and Width of scattering region (the superconductor) L = 20 W = 10 def make_system(): lat = kwant.lattice.general([(sqrt(3)*0.5,0.5),(0,1)],[(0,0),(1/(2*sqrt(3)),1/2)],norbs = 4) a, b = lat.sublattices def channel(pos): x, y = pos return -W <= y <= W and -L <= x <= L syst = kwant.Builder() syst[lat.shape(channel, (0, 0))] = Ef_S * -1 * np.kron(tau_z, s_0) syst[lat.neighbors()] = t * -1 * np.kron(tau_z, s_0) syst.eradicate_dangling() def lead0_shape(pos): # left lead x, y = pos return -W <=y <=W sym_up = kwant.TranslationalSymmetry(lat.vec((-2,1))) # Symmetry along the x-axis # Normal graphene lead-0 lead0 = kwant.Builder(sym_up, conservation_law=np.kron(tau_z, s_0), particle_hole=np.kron(tau_x, s_0)) lead0[lat.shape(lead0_shape, (0, 0))] = Ef_S * -1 * np.kron(tau_z, s_0) lead0[lat.neighbors()] = t * -1 * np.kron(tau_z, s_0) ######################################################################################################################## def lead1_shape(pos): # right x, y = pos return -W <=y <=W sym_down = kwant.TranslationalSymmetry(lat.vec((2, -1))) # Superconducting graphene lead-1 lead1 = kwant.Builder(sym_down) lead1[lat.shape(lead1_shape, (0, 0))] = (Ef_S + V0) * -1 * np.kron(tau_z, s_0) + Delta * np.kron(tau_z, s_0) lead1[lat.neighbors()] = t * -1 * np.kron(tau_z, s_0) return syst, [lead0, lead1] def plot_conductance(syst, energies): # Compute conductance data = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) # Conductance is N - R_ee + R_he data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0))) plt.figure() plt.plot(energies, data,'*-r') plt.xlabel("energy [t]") plt.ylabel("conductance [e^2/h]") plt.show() #Plot the band in the lead #Make my system and attach the leads def main(): system, leads = make_system() # plot_system(system, a) # Attach the leads to the system. for lead in leads: system.attach_lead(lead) # # Plot system with leads # def family_colors(site): # return 0 if site.family == a else 1 # # kwant.plot(system, site_color=family_colors, site_lw=0.1, lead_site_lw=0, colorbar=False) # Finalize the system. kwant.plot(system,fig_size=(5,5)) system = system.finalized() # #Evaluate the band in the lead L_1 = leads[0].finalized() # bands = kwant.physics.Bands(L_1) # #Plot the band in the lead # ks = np.linspace(-pi,pi,num=101) # energies=[bands(k) for k in ks] # f, ax = plt.subplots() # ax.plot(ks,energies) # plt.show() plot_conductance(system, energies=[0.01 * i+0.00001 for i in range(200)]) main()

1 0

Methods for calculating Josephson Current ?
by guilhermewells＠gmail.com 06 Jun '23

06 Jun '23

Hi, I'm new to Kwant and I've been looking for calculating Josephson Current with it. The past three weeks I've been reading a lot the tutorials, the mailing list and some online courses. The problem is that I can't get a straight answer on where it is possible to calculate Josephson current reallistically. The most complete answer was given by J. Weston here [1], the problem is that the methods suggested are also not so clear, since I coundn't find tutorials using properly Self-Energy and Green's function with Kwant. Most of the topics are questions, so they usually are of people with code not working, so it doesn't point out the proper way to use. So, I'd like to know if I'm missing some tutorials with these subjects, or some online courses to understand this. The documentation itself is not enough for the applications (at least for me). I'm learning kwant alone, there is no one around in my department that even have heard of it. Right now I'm got a discontinuous current (only at the junction interface) through my junction using directly "kwant.operator.Current" and "plotter", which I think is to be expected due to the issue of the Andreev bound states. Thanks, I'm looking foward to use kwant at higher level [1] https://mail.python.org/archives/list/kwant-discuss@python.org/thread/6IP2Q… -- Guilherme Vieira Phd Student at Centro Brasileiro de Pesquisas Físicas Brazil, RJ

2 1

modelling chiral channels
by Sergey Slizovskiy 05 Jun '23

05 Jun '23

Hi, I am thinking if it's possible to use Kwant to compute the conductance matrix for a network of chiral 1D channels? Looks very similar to what kwant normally does, I suppose, but kwant does not allow for non-symmetric hoppings. Any ideas? Thank you, Sergey

2 1