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
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
-- Abbout Adel
participants (2)
-
Abbout Adel
-
mondal18@iitg.ac.in