Dear kwant developers and kwant community,
I have been using kwant for the past three months, and I have
encountered a strange outcome. It may be just a mistake from my
part, but after some testing the outcome makes little sense to me. I
would really appreciate your help and/or input.
Using the script below, I have computed the wavefunction of a
(quite) wide graphene armchair lead with a huge (50 t, ~"infinite")
mass near the boundary. It turns out that the outcome strongly
depends on the x-coordinates of the system. With the variable xShift
set to zero, the script produces a ribbon with x-coordinates between
-W/2 and W/2. Setting xShift to a positive value shifts this to the
left.
Now, it turns out that setting the variable 'xShift' to zero
produces a wavefunction that differs quite a lot from the
wavefunction produced by setting xShift to other values, e.g.
4*sqrt(3). Furthermore, the outcomes are different for kwant
versions 1.0.0 and 1.0.1. Going through the git log, it was not
obvious to me which change could cause these different outcomes. I
tested one (much bigger) system on kwant 1.0.3, and it seemed to
give the same outcome as 1.0.1
Even stranger, running the same system with xShift=1.5*sqrt(3) on
kwant version 1.0.1 sometimes gives two different outcomes, although
this does not seem to be very reproducable. Could this mean that
there is some instability in the computation? (There are no error
messages shown)
This effect does not seem to occur (as far as I have checked) for
- Armchair ribbons without a mass term (i.e. setting mAmpEdge to 0
in my script)
- Zigzag ribbons with a mass term (I can send the code if that would
help)
- Zigzag ribbons without a mass term
Could this mean that this strange outcome is due to the peculiarity
of this system?
I have run kwant.test() for both kwant versions (albeit on two
different systems), which gives two errors for the system running
kwant version 1.0.0, see the output below, and no errors for the
system running kwant version 1.0.1.
Thank you in advance for your input.
Best regards,
Koen Reijnders
PhD-student
Radboud University
Nijmegen, The Netherlands
=================== Script ======================================
import kwant
import numpy
import tinyarray
from matplotlib import pyplot
from math import pi, sqrt, exp
import math
lat = kwant.lattice.general([(sqrt(3), 0), (sqrt(3)/2., 3/2.)],
[(0, 0), (0, 1)])
a, b = lat.sublattices
W=(6*100-1)*0.5*sqrt(3); # width of the lead, in units of aCC
# Of the form (6*N-1)*0.5*sqrt(3), N INTEGER, to make the lead
metallic
xShift=0*sqrt(3) # shift in the position of the lead
En=0.198/3. # energy, in units of t
incomingModeIdx=0; # index of the mode in the lead
noModes=11; # number of modes to be printed
mAmpEdge=50; # mass introduced at the edge
ledge=30; # length over which the mass is introduced at the edge
# Creates the lead
def make_lead(left_edge,right_edge):
sym = kwant.TranslationalSymmetry((0, -3))
lead = kwant.Builder(sym)
def lead_shape(pos):
x, y = pos
return left_edge < x < right_edge
# Infinite mass potential
def massEdge(site):
(x, y) = site.pos
msign = 1 if site.family == a else -1
medge =
msign*mAmpEdge*(2-numpy.tanh((x-left_edge)/ledge)+numpy.tanh((x-right_edge)/ledge))
return medge
lead[lat.shape(lead_shape, (1,1))] = massEdge
lead[lat.neighbors()] = -1.0
return lead
# Determines the spectrum of the lead. Time-consuming step
def determine_lead_spectrum(alead):
finlead=alead.finalized()
leadSpectrum=finlead.modes(En)[0]
return finlead, leadSpectrum
# Orders the wavefunction in the lead. Returns the wavefunction
(ordered)
def extract_lead_wavefunctions(finlead, propModes,
imode=incomingModeIdx):
sizeMomenta=numpy.size(propModes.momenta)
modeIdx=sizeMomenta/2+imode
### THE FOLLOWING CODE ONLY WORKS FOR A SINGLE LEAD CELL!!! ###
# Taken from the Kwant source code: kwant.plotter.sys_leads_pos
and related. #
num_lead_cells=1
sites=[(site, i) for site in xrange(finlead.cell_size) for i in
xrange(num_lead_cells)]
pos = numpy.array(tinyarray.array([finlead.pos(i[0]) for i in
sites]))
matSitesWF=numpy.zeros((len(pos),4))
matSitesWF[:,:2]=pos
matSitesWF[:,2]=numpy.real(propModes.wave_functions[:,modeIdx])
matSitesWF[:,3]=numpy.imag(propModes.wave_functions[:,modeIdx])
minpos=numpy.min(pos[:,0])
Xstart=math.floor(minpos/sqrt(3))*sqrt(3)
tol=1e-8
safetyMargin=2;
matSitesOrdered=numpy.zeros((len(pos)/4+safetyMargin,4,4))
for i in xrange(len(pos)):
rowI=math.floor((matSitesWF[i,0]-Xstart+tol)/sqrt(3))
if matSitesWF[i,1]==-1.5:
indJ=0;
elif matSitesWF[i,1]==-0.5:
indJ=1
elif matSitesWF[i,1]==0:
indJ=2
elif matSitesWF[i,1]==1.0:
indJ=3
else:
print "Something wrong"
matSitesOrdered[rowI,indJ,:]=matSitesWF[i,:]
return matSitesOrdered
# Write to file
def
write_wavefunctions(matSitesOrdered,width,energy,incomingModeIdx):
filenameStr="LeadWF_InfMassAC_W" + str(width) + "_En" +
str(energy) + "_mode" + str(incomingModeIdx) + ".dat"
with file(filenameStr,"w") as outfile:
for slice2d in matSitesOrdered:
numpy.savetxt(outfile,slice2d)
# Load from file
def load_wfdata(width,energy,incomingModeIdx):
filenameStr="LeadWF_InfMassAC_W" + str(width) + "_En" +
str(energy) + "_mode" + str(incomingModeIdx) + ".dat"
return numpy.loadtxt(filenameStr)
# Plot functions for the lead data
def plot_abs(WFData):
sitesAback = WFData[0::4,:]
sitesBback = WFData[1::4,:]
sitesAfront = WFData[2::4,:]
sitesBfront = WFData[3::4,:]
fig=pyplot.figure()
ax=pyplot.subplot(111)
line1=ax.plot(sitesAback[:,0],numpy.sqrt(numpy.square(sitesAback[:,2])+numpy.square(sitesAback[:,3])),'ro',label="A
back")
line2=ax.plot(sitesAfront[:,0],numpy.sqrt(numpy.square(sitesAfront[:,2])+numpy.square(sitesAfront[:,3])),'r^',label="A
front")
line3=ax.plot(sitesBback[:,0],numpy.sqrt(numpy.square(sitesBback[:,2])+numpy.square(sitesBback[:,3])),'bo',label="B
back")
line4=ax.plot(sitesBfront[:,0],numpy.sqrt(numpy.square(sitesBfront[:,2])+numpy.square(sitesBfront[:,3])),'b^',label="B
front")
ax.set_ylabel("|Psi|")
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1, box.width,
box.height * 0.9])
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.05),
fancybox=True, shadow=True, ncol=4)
pyplot.show()
# Main method
def main():
# create the lead and determine its spectrum
lead=make_lead(-W/2.0-xShift,W/2.0-xShift)
flead,leadSpec=determine_lead_spectrum(lead)
#Optional: writing modes to file for later use
# for modePrint in xrange(noModes):
# wflead=extract_lead_wavefunctions(flead,leadSpec,modePrint)
# write_wavefunctions(wflead,W,En,modePrint)
# write_momenta(leadSpec,W,En)
wflead=extract_lead_wavefunctions(flead,leadSpec,incomingModeIdx)
wfleadflat=numpy.array([item for sublist in wflead for item in
sublist])
plot_abs(wfleadflat)
if __name__ == "__main__":
main()
===========================================================
System 1: Ubuntu Linux 14.04
Python 2.7.6 (default, Jun 22 2015, 17:58:13)
[GCC 4.8.2] on linux2
>>> import kwant
>>> kwant.version.version
'1.0.1'
>>> kwant.test()
........................SSS................SSSSSSSSSSSS.........................................
----------------------------------------------------------------------
Ran 96 tests in 8.929s
OK (SKIP=15)
True
System 2: Ubuntu Linux 12.04
Python 2.7.3 (default, Jun 22 2015, 19:33:41)
[GCC 4.6.3] on linux2
>>> import kwant
>>> kwant.version.version
'1.0.0'
>>> kwant.test()
.......................................................E.......................E....
======================================================================
ERROR: Failure: AttributeError ('module' object has no attribute
'umfpack')
----------------------------------------------------------------------
Traceback (most recent call last):
File "/usr/lib/python2.7/dist-packages/nose/loader.py", line 390,
in loadTestsFromName
addr.filename, addr.module)
File "/usr/lib/python2.7/dist-packages/nose/importer.py", line 39,
in importFromPath
return self.importFromDir(dir_path, fqname)
File "/usr/lib/python2.7/dist-packages/nose/importer.py", line 86,
in importFromDir
mod = load_module(part_fqname, fh, filename, desc)
File
"/usr/lib/python2.7/dist-packages/kwant/solvers/tests/test_sparse.py",
line 10, in <module>
from kwant.solvers.sparse import smatrix, greens_function,
ldos, wave_function
File "/usr/lib/python2.7/dist-packages/kwant/solvers/sparse.py",
line 18, in <module>
umfpack = linsolve.umfpack
AttributeError: 'module' object has no attribute 'umfpack'
======================================================================
ERROR: kwant.tests.test_plotter.test_plot
----------------------------------------------------------------------
Traceback (most recent call last):
File "/usr/lib/python2.7/dist-packages/nose/case.py", line 197, in
runTest
self.test(*self.arg)
File
"/usr/lib/python2.7/dist-packages/kwant/tests/test_plotter.py", line
102, in test_plot
fig = plot(sys, site_color=color, cmap='binary', file=out)
File "/usr/lib/python2.7/dist-packages/kwant/plotter.py", line
1361, in plot
return output_fig(fig, file=file, show=show)
File "/usr/lib/python2.7/dist-packages/kwant/plotter.py", line
663, in output_fig
canvas.print_figure(file, *savefile_opts[0], **savefile_opts[1])
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/backend_bases.py",
line 2192, in print_figure
**kwargs)
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.py",
line 513, in print_png
FigureCanvasAgg.draw(self)
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.py",
line 461, in draw
self.figure.draw(self.renderer)
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/artist.py", line
59, in draw_wrapper
draw(artist, renderer, *args, **kwargs)
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/figure.py", line
1079, in draw
func(*args)
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/artist.py", line
59, in draw_wrapper
draw(artist, renderer, *args, **kwargs)
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/axes/_base.py",
line 2092, in draw
a.draw(renderer)
File "/usr/lib/python2.7/dist-packages/kwant/plotter.py", line
117, in draw
return collections.Collection.draw(self, renderer)
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/artist.py", line
59, in draw_wrapper
draw(artist, renderer, *args, **kwargs)
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/collections.py",
line 301, in draw
self._offset_position)
File
"/usr/local/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.py",
line 119, in draw_path_collection
return self._renderer.draw_path_collection(*kl, **kw)
ValueError: Transforms must be a Nx3x3 numpy array
----------------------------------------------------------------------
Ran 84 tests in 18.351s
FAILED (errors=2)
False