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
