Dear all,
I am new to kwant. I am trying to calculate the current density for different
widths of zigzag graphene nanoribbon. But I am getting strange results for
widths above 22 angstroms, there is slightly back propagation of current. I
want to know, is there any mistake in my code or anything else to be considered
for the calculation of current density. Here i am attaching the code with
current density result for width 22 and 24 angstroms. In 24 angstroms, there is
back propagation of current. Please help me.
--------Code ----
import kwant
import numpy as np
from matplotlib import pyplot
import matplotlib.pyplot as plt
import numpy.linalg as npl
import scipy.sparse.linalg as sla
import os
pi=np.pi
############################################## structure
########################
Wnr=22
lnr=180
############# empty system kwant
lat=kwant.lattice.honeycomb(a=2.46,norbs=1)
a,b=lat.sublattices
########### functions for structure build
def onsite(site, voltage):
x, y = site.pos
return voltage
def rect(pos):
x,y=pos
return 0<x<=lnr and 0<y<=Wnr
################################ main structure i.e. rectangle
####################
model1 = kwant.Builder()
model1[lat.shape(rect,(1,1))] = 0
model1[lat.neighbors()] =2.64
model1.eradicate_dangling()
kwant.plot(model1)
############### functions for lead structure build
################################
def lead1_shape(pos):
x, y = pos
return 0<y <=Wnr
############first lead codes##########
sym = kwant.TranslationalSymmetry(lat.vec((-1,0))) #from the left
sym.add_site_family(lat.sublattices[0], other_vectors=[(-1, 2)])
sym.add_site_family(lat.sublattices[1], other_vectors=[(-1, 2)])
lead = kwant.Builder(sym)
lead[lat.shape(lead1_shape, (0,1))]=onsite
lead[lat.neighbors()] =2.64
#########attaching the lead to model1
model1.attach_lead(lead)
model1.attach_lead(lead.reversed()) #second lead attach in reverse
manner
model1=model1.finalized()
kwant.plot(model1)
####################### calculations ############################
###functions for calculations
params = dict(voltage=0)
kwant.plotter.bands(model1.leads[1],params=params) #bandstructure of
lead
def plot_conductance(sys, energies,params):
# Compute transmission as a function of energy
data = []
for energy in energies:
smatrix = kwant.smatrix(sys, energy, params=params)
data.append(smatrix.transmission(0, 1))
pyplot.figure()
pyplot.plot(energies, data)
#pyplot.xlim([-2,2])
#pyplot.ylim([0,6])
pyplot.xlabel("energy [t]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()
return data
##########DOS KPM#################
rho = kwant.kpm.SpectralDensity(model1,params=params)
energies, densities = rho.energies, rho.densities
plt.figure()
plt.plot(energies, densities)
#####################################
#######conductance w.r.t energy##############
data=plot_conductance(model1,energies, params)
#######################
###########current density########
J_0 = kwant.operator.Current(model1)
wf = kwant.wave_function(model1,energy=1,params=params)
wfs=wf(0)
wfs_of_lead_2 = wf(1)
psi = wf(0)[0]
current = J_0(psi)
kwant.plotter.current(model1, current, colorbar=True)
################################################
please find the current density result in the link
https://drive.google.com/file/d/1j8A6yZbZwaR1yAi3gKnuZalJYTJ8FQHc/view?usp=sharing
https://drive.google.com/file/d/155-PxKLLCtGewqQL2fuitGWPI6KCReik/view?usp=sharing
