Dear Koen, Sorry for the delayed reply. The strange behavior that you observe isn't related to Kwant, and instead it's a property of graphene boundary conditions. Basically when people use the term "infinite mass" they mean mass much larger than Fermi energy but still much smaller than the hopping strength. For more details about boundary conditions in graphene you can read for example this paper: http://arxiv.org/abs/0710.2723 The run-dependent modes correspond to degenerate solutions in which case the basis choice becomes arbitrary. In that case physical observables won't change regardless of the basis choice. Finally the test failures don't site anything bad. Best, Anton On Tue, Jul 14, 2015, 14:41 Koen Reijnders <K.Reijnders@science.ru.nl> wrote:
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