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