Dear developers, I was trying to get the plot of currents in a graphene ribbon. I am getting the error "Number of orbitals not defined", which is done by specifying the norbs. However, in my code, I am unable to specify the norbs. Another thing is that the interface between the right lead (and left leads) and the system is not parallel to the y axis. What are the possible ways to solve such problems. The code I am using is the following. #################################################### from math import sqrt import random import itertools as it import tinyarray as ta import numpy as np import matplotlib.pyplot as plt import kwant class Honeycomb(kwant.lattice.Polyatomic): """Honeycomb lattice with methods for dealing with hexagons""" def __init__(self, name=''): prim_vecs = [[0.5, sqrt(3)/2], [1, 0]] # bravais lattice vectors # offset the lattice so that it is symmetric around x and y axes basis_vecs = [[-0.5, -1/sqrt(12)], [-0.5, 1/sqrt(12)]] super(Honeycomb, self).__init__(prim_vecs, basis_vecs, name) self.a, self.b = self.sublattices def hexagon(self, tag): """ Get sites belonging to hexagon with the given tag. Returns sites in counter-clockwise order starting from the lower-left site. """ tag = ta.array(tag) # a-sites b-sites deltas = [(0, 0), (0, 0), (1, 0), (0, 1), (0, 1), (-1, 1)] lats = it.cycle(self.sublattices) return (lat(*(tag + delta)) for lat, delta in zip(lats, deltas)) def hexagon_neighbors(self, tag, n=1): """ Get n'th nearest neighbor hoppings within the hexagon with the given tag. """ hex_sites = list(self.hexagon(tag)) return ((hex_sites[(i+n)%6], hex_sites[i%6]) for i in range(6)) def random_placement(builder, lattice, density): """ Randomly selects hexagon tags where adatoms can be placed. This avoids the edge case where adatoms would otherwise be placed on incomplete hexagons at the system boundaries. """ assert 0 <= density <= 1 system_sites = builder.sites() def hexagon_in_system(tag): return all(site in system_sites for site in lattice.hexagon(tag)) # get set of tags for sites in system (i.e. tags from only one sublattice) system_tags = (s.tag for s in system_sites if s.family == lattice.a) # only allow tags which have complete hexagons in the system valid_tags = list(filter(hexagon_in_system, system_tags)) N = int(density * len(valid_tags)) total_hexagons=len(valid_tags) valuef=random.sample(valid_tags, N) return valuef def ribbon(W, L): def shape(pos): return (-L <= pos[0] <= L and -W <= pos[1] <= W) return shape ## Pauli matrices ## sigma_0 = ta.array([[1, 0], [0, 1]]) # identity sigma_x = ta.array([[0, 1], [1, 0]]) sigma_y = ta.array([[0, -1j], [1j, 0]]) sigma_z = ta.array([[1, 0], [0, -1]]) ## Hamiltonian value functions ## def onsite_potential(site, params): return params['ep'] * sigma_0 def potential_shift(site, params): return params['mu'] * sigma_0 def kinetic(site_i, site_j, params): return -params['gamma'] * sigma_0 def rashba(site_i, site_j, params): d_ij = site_j.pos - site_i.pos rashba = 1j * params['V_R'] * (sigma_x * d_ij[1] - sigma_y * d_ij[0]) return rashba + kinetic(site_i, site_j, params) def spin_orbit(site_i, site_j, params): so = 1j * params['V_I'] * sigma_z return so lat = Honeycomb() A, B = lat.sublattices pv1, pv2 = lat.prim_vecs ysym = kwant.TranslationalSymmetry(pv2 - 2*pv1) # lattice symmetry in -y direction xsym = kwant.TranslationalSymmetry(-pv2) # lattice symmetry in -x direction # adatom lattice, for visualization only adatom_lat = kwant.lattice.Monatomic(lat.prim_vecs, name='adatom') def site_size(s): return 0.25 if s.family in lat.sublattices else 0.3 def site_color(s): if s.family==lat.sublattices[0]: return 'w' elif s.family==lat.sublattices[1]: return 'k' else: return 'g' def create_lead_h(W, sym1, axis=(0, 0)): lead = kwant.Builder(sym1) lead[lat.wire(axis, W)] = 0. * sigma_0 lead[lat.neighbors(1)] = kinetic return lead def create_ribbon(W, L, adatom_density=0.2, show_adatoms=False): ## scattering region ## sys = kwant.Builder() sys[lat.shape(ribbon(W, L), (0, 0))] = onsite_potential sys[lat.neighbors(1)] = kinetic adatoms = random_placement(sys, lat, adatom_density) sys[(lat.hexagon(a) for a in adatoms)] = potential_shift sys[(lat.hexagon_neighbors(a, 1) for a in adatoms)] = rashba sys[(lat.hexagon_neighbors(a, 2) for a in adatoms)] = spin_orbit if show_adatoms: # no hoppings are added so these won't affect the dynamics; purely cosmetic sys[(adatom_lat(*a) for a in adatoms)] = 0 ## leads ## leads = [create_lead_h(W, xsym)] leads += [lead.reversed() for lead in leads] # right lead for lead in leads: sys.attach_lead(lead) return sys def plot_currents(syst, currents): fig, axes = plt.subplots(1, len(currents)) if not hasattr(axes, '__len__'): axes = (axes,) for ax, (title, current) in zip(axes, currents): kwant.plotter.current(syst, current, ax=ax, colorbar=False) ax.set_title(title) plt.show() def currents(sys, params): syst = sys.finalized() wf = kwant.wave_function(syst, energy=-1, args=(params,)) psi = wf(0)[0] # even (odd) indices correspond to spin up (down) up, down = psi[::2], psi[1::2] print(psi.shape) density = np.abs(up)**2 + np.abs(down)**2 # spin down components have a minus sign spin_z = np.abs(up)**2 - np.abs(down)**2 # spin down components have a minus sign spin_y = 1j * (down.conjugate() * up - up.conjugate() * down) rho = kwant.operator.Density(syst) rho_sz = kwant.operator.Density(syst, sigma_z) rho_sy = kwant.operator.Density(syst, sigma_y) # calculate the expectation values of the operators with 'psi' density = rho(psi) spin_z = rho_sz(psi) spin_y = rho_sy(psi) plot_densities(syst, [ ('$σ_0$', density), ('$σ_z$', spin_z), ('$σ_y$', spin_y), ]) J_0 = kwant.operator.Current(syst) J_z = kwant.operator.Current(syst, sigma_z) J_y = kwant.operator.Current(syst, sigma_y) # calculate the expectation values of the operators with 'psi' current = J_0(psi) spin_z_current = J_z(psi) spin_y_current = J_y(psi) plot_currents(syst, [ ('$J_{σ_0}$', current), ('$J_{σ_z}$', spin_z_current), ('$J_{σ_y}$', spin_y_current), ]) if __name__ == '__main__': params = dict(gamma=1., ep=0, mu=0., V_I=0.007, V_R=0.0165) # In adatoms W=11 L=32 adatom_density=0.1 sys = create_ribbon(W, L, adatom_density, show_adatoms=True) # fig = plt.figure() ax = plt.subplot() ax = kwant.plot(sys, site_color=site_color, site_lw=0.05, site_size=site_size,lead_site_lw=0, dpi = 150) # ax.savefig('system2.png', bbox_inches = 'tight', dpi = 200) currents(sys, params) ################################################### I shall look forward to your kind reply. Thanking you. Sincerely yours, Sayan Mondal, Research scholar, IIT Guwahati, India