Dear Koustav, Your example is puzzling me and I could not find the mistake! If you change the spin orbit coupling lambda_so = 0.1 to higher values like 0.6 and above it works! If you change the value of the width to 30 sqrt(1/3) or lower it works!
Please post the solution if you find it. Adel On Mon, Jun 14, 2021 at 6:24 AM <170070...@iitb.ac.in> wrote: > Hi, > The code I am using is pasted below. This gives the above mentioned > RunTime Error. > > import kwant > from math import pi,sqrt > import numpy as np > from matplotlib import pyplot > from kwant.digest import uniform > from types import SimpleNamespace > > pauli_z = np.array([[1,0],[0,-1]]) > sin_30, cos_30 = (1 / 2, sqrt(3) / 2) > graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)], > [(0, 0), > (0, 1 / sqrt(3))]) > > a, b = graphene.sublattices > > def lin0(y,W,jw) : > if y < -jw : > return -2 > elif -jw <= y < jw : > return 2*y/jw > else : > return 2 > > def lin1(y,W,jw) : > if y < -W/6-jw : > return -2 > elif -W/6-jw <= y < -W/6+jw : > return (y+W/6)/jw - 1 > elif -W/6+jw <= y < W/6-jw : > return 0 > elif W/6-jw <= y < W/6+jw : > return (y-W/6)/jw + 1 > else : > return 2 > > def lin2(y,W,jw) : > if y < -jw : > return 0 > elif -jw <= y < jw : > return y/jw+1 > else : > return 2 > > > def make_system(W = 10*sqrt(3), L = 10, delta = 0, t = 1.6, lambda_so = 0) > : > > lambda_so = 1j*lambda_so/(3*sqrt(3)) > t_nn1_a = 1 * lambda_so * pauli_z # [ 1, 0] > t_nn1_b = -1 * lambda_so * pauli_z > t_nn2_a = -1 * lambda_so * pauli_z # [ 0, 1] > t_nn2_b = 1 * lambda_so * pauli_z > t_nn3_a = -1 * lambda_so * pauli_z # [ 1, -1] > t_nn3_b = 1 * lambda_so * pauli_z > > def channel(pos): > x, y = pos > return (0 <= x <= L) and (-W/2 < y <= W/2) > > syst = kwant.Builder() > > del_fn = lambda y,W : lin1(y,W,W/20) > > def potential(site, U, U_disorder, Mex): > (x, y) = site.pos > salt = 0 > d = -1 > if (site.family == a) : > d = 1 > term1 = d*delta*del_fn(y,W)*np.eye(2) > term2 = U*np.eye(2) > term3 = Mex*pauli_z > term4 = U_disorder * (uniform(repr(site), repr(salt)) - > 0.5) * np.eye(2) > return term1 + term2 + term3 + term4 > > > def dummy(site, Mex): > (x, y) = site.pos > d = -1 > if (site.family == a) : > d = 1 > term1 = d*delta*del_fn(y,W)*np.eye(2) > term2 = Mex*pauli_z > return term1 + term2 > > > syst[graphene.shape(channel, (0, 0))] = potential > hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) > syst[[kwant.builder.HoppingKind(*hopping) for hopping in > hoppings]] = -t*np.eye(2) > syst[kwant.builder.HoppingKind((1, 0), a, a)] = t_nn1_a > syst[kwant.builder.HoppingKind((1, 0), b, b)] = t_nn1_b > syst[kwant.builder.HoppingKind((0, 1), a, a)] = t_nn2_a > syst[kwant.builder.HoppingKind((0, 1), b, b)] = t_nn2_b > syst[kwant.builder.HoppingKind((1, -1), a, a)] = t_nn3_a > syst[kwant.builder.HoppingKind((1, -1), b, b)] = t_nn3_b > > # left lead > sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0))) > > def lead0_shape(pos): > x, y = pos > return (-W/2 < y <= W/2) > > lead0 = kwant.Builder(sym0) > lead0[graphene.shape(lead0_shape, (0, 0))] = np.zeros((2,2)) > lead0[[kwant.builder.HoppingKind(*hopping) for hopping in > hoppings]] = -t*np.eye(2) > lead0[kwant.builder.HoppingKind((1, 0), a, a)] = t_nn1_a > lead0[kwant.builder.HoppingKind((1, 0), b, b)] = t_nn1_b > lead0[kwant.builder.HoppingKind((0, 1), a, a)] = t_nn2_a > lead0[kwant.builder.HoppingKind((0, 1), b, b)] = t_nn2_b > lead0[kwant.builder.HoppingKind((1, -1), a, a)] = t_nn3_a > lead0[kwant.builder.HoppingKind((1, -1), b, b)] = t_nn3_b > > # right lead > sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0))) > > def lead1_shape(pos): > x, y = pos > return (-W/2 < y <= W/2) > > lead1 = kwant.Builder(sym1) > lead1[graphene.shape(lead1_shape, (0, 0))] = np.zeros((2,2)) > lead1[[kwant.builder.HoppingKind(*hopping) for hopping in > hoppings]] = -t*np.eye(2) > lead1[kwant.builder.HoppingKind((1, 0), a, a)] = t_nn1_a > lead1[kwant.builder.HoppingKind((1, 0), b, b)] = t_nn1_b > lead1[kwant.builder.HoppingKind((0, 1), a, a)] = t_nn2_a > lead1[kwant.builder.HoppingKind((0, 1), b, b)] = t_nn2_b > lead1[kwant.builder.HoppingKind((1, -1), a, a)] = t_nn3_a > lead1[kwant.builder.HoppingKind((1, -1), b, b)] = t_nn3_b > > # dummy lead > dum_sym = kwant.TranslationalSymmetry(graphene.vec((1, 0))) > > def dum_lead_shape(pos): > x, y = pos > return (-W/2 < y <= W/2) > > dum_lead = kwant.Builder(dum_sym) > dum_lead[graphene.shape(dum_lead_shape, (0, 0))] = dummy > dum_lead[[kwant.builder.HoppingKind(*hopping) for hopping in > hoppings]] = -t*np.eye(2) > dum_lead[kwant.builder.HoppingKind((1, 0), a, a)] = t_nn1_a > dum_lead[kwant.builder.HoppingKind((1, 0), b, b)] = t_nn1_b > dum_lead[kwant.builder.HoppingKind((0, 1), a, a)] = t_nn2_a > dum_lead[kwant.builder.HoppingKind((0, 1), b, b)] = t_nn2_b > dum_lead[kwant.builder.HoppingKind((1, -1), a, a)] = t_nn3_a > dum_lead[kwant.builder.HoppingKind((1, -1), b, b)] = t_nn3_b > > return syst, [lead0, lead1], dum_lead > > W = 40*sqrt(3) > L = 40 > t = 1.3 > E = 0.5 > U = E > lambda_so = 0.1 > delta = -lambda_so > U_disorder = 4*delta > params = dict(U = U, U_disorder = U_disorder, Mex = 0) > > > syst, leads, dum_lead = make_system(W = W, L = L, delta = delta, t = t, > lambda_so = lambda_so) > > > for lead in leads: > syst.attach_lead(lead) > > syst = syst.finalized() > > local_dos = kwant.ldos(syst,energy=E,params=params) > kwant.plotter.map(syst, local_dos[1::2], num_lead_cells=0, a=1/sqrt(3)) > -- Abbout Adel
