Problems in getting the current
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
Dear Sayan, You have to specify the number of orbitals in the lattice: norbs=2. You also have to make sure that the number of sites goes well with the shape of the wavefunction. Avoid this sentence that changes the number of sites: if show_adatoms: # no hoppings are added so these won't affect the dynamics; purely cosmetic pass #sys[(adatom_lat(*a) for a in adatoms)] = 0 I strongly suggest using params instead of args. The way you used params is not the standard way. You forgot to use the params in the operator current and density. With these corrections, the following code works. I hope this helps, Adel ################################################################# import kwant import tinyarray as ta from numpy import sqrt import itertools as it from random import random,sample import numpy as np import random from matplotlib import pyplot as plt class Honeycomb(kwant.lattice.Polyatomic): """Honeycomb lattice with methods for dealing with hexagons""" def __init__(self, name='', norbs=2): 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,norbs=2) 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=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',norbs=2) 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 pass#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 plot_densities(syst, densities): fig, axes = plt.subplots(1, len(densities), figsize=(13, 10)) for ax, (title, rho) in zip(axes, densities): kwant.plotter.density(syst, rho, ax=ax) 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) print(psi) rho = kwant.operator.Density(syst,sigma_0) 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' print ("shape= ",psi.shape) print("sites=:", len(list(syst.sites)),"double ",2*len(list(syst.sites))) # psi=psi[::2] density = rho(psi,args=(params,)) spin_z = rho_sz(psi,args=(params,)) spin_y = rho_sy(psi,args=(params,)) 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,args=(params,)) spin_z_current = J_z(psi,args=(params,)) spin_y_current = J_y(psi,args=(params,)) 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) On Mon, May 9, 2022 at 12:53 PM <mondal18@iitg.ac.in> wrote:
-- Abbout Adel
participants (2)
-
Abbout Adel
-
mondal18@iitg.ac.in