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
