Use better heuristic for colorscale limits
Hello! I am here to report a very unexpected behavior I experience when plotting wavefunctions of graphene pn-junctions.
The curious thing that happens is that when plotting the wavefunction incoming from the second graphene valley (more on that below), the plotting complete fails. What is shown instead is a uniform blob of color, normally the lowest value of the colormap.
This cannot be physically correct, as it shows that the wavefunction simply doesn't exist. It is extremely curious that this happens only on one of the two valleys. Graphene zigzag nanoribbons have wavefunctions localized on one of two possible (near)-Dirac valleys in momentum space. The wavefunctions from one valley are plotted properly but from the other not.
Here is my example, where I plot mode 0
and mode 8
, which are the lowest modes from each valley.
It gets even weirder when you realize that the current
plotter works perfectly fine!
Below you can find the (quite long) source code so that you can immediately replicate my findings.
I took the liberty to post this as an issue and not ask about it in the mailing list, because I sincerely believe it is a code related issue and not a fault at my understanding. Please forgive me if this is not the case.
Notice that I have indeed used these wavefunctions for research purposes, and everything seems to be "ok", in the sense that my research shows reasonable results using the wavefunctions from both valleys, making me even further suspect that this is a plotting issue (somehow?).
import kwant
import numpy as np
from numpy import sqrt, exp, pi, cos, sin
import matplotlib.pyplot as plt
from matplotlib.pyplot import *
from mpl_toolkits.mplot3d import Axes3D
import scipy
import types
DEFAULT_PAR = types.SimpleNamespace(B = 0)
t = 2.8 #hopping in eV
a0 = 0.142*sqrt(3) #graphene lattice constant
e_charge = 1.6021766208e-19 # SI units
hbar = 1.054571800e-34 # plancks constant (SI Units)
prefac = e_charge/hbar
# hbar = 6.5821195141e-16 #planck in eV*s
vF = 1e6 * 1e9 #vF in nm/s
hbarvF = 0.6582119514 #in units of eV*nm
def make_graphene(shape, pot=0.0, init = (0,0), a0 = a0, t=t):
glat = kwant.lattice.honeycomb(a0, norbs=1)
alat, blat = glat.sublattices
gsys = kwant.Builder()
def hopping(sitei, sitej, p):
# FROM sitej TO sitei
xi, yi = sitei.pos
xj, yj = sitej.pos
return -t * exp(1j * prefac * p.B * (xi - xj) * (yi + yj) * 0.5*1e-18)
gsys[glat.shape(shape, init)] = pot # make the actual graphene
gsys[glat.neighbors()] = hopping
gsys.eradicate_dangling()
return gsys, glat, alat, blat
def make_graphene_pnjunction(length, wid, V0=0.1, pottype="linear", w=2.0):
V02 = V0/2
def shape(pos): # shape of junction
x, y = pos
return 0 <= x <= length and 0 <= y <= wid
def spacepot_linear(x, y):
if x <= length/2 - w/2:
return -V02
elif x >= length/2 + w/2:
return V02
else:
return (x-length/2)*V0/w
def spacepot_tanh(x, y):
return np.tanh((x-length/2)/(w/2))*V0/2
def spacepot_exp(x, y):
return V0/(1 + np.exp(-(x-length/2)/(w/2))) - V02
if pottype == "linear":
def potential(site, p): #potential must use parameters to allow B
x, y = site.pos
return spacepot_linear(x,y)
elif pottype == "tanh":
def potential(site, p):
x, y = site.pos
return spacepot_tanh(x,y)
elif pottype == "exp":
def potential(site, p):
x, y = site.pos
return spacepot_exp(x,y)
gsys, glat, alat, blat = make_graphene(shape, pot = potential)
leftlead = make_zigzag_lead(wid, pot = -V02)
rightlead = make_zigzag_lead(wid, init=(length,0), pot= V02, direction=+1)
gsys.attach_lead(leftlead);
gsys.attach_lead(rightlead);
fsys = gsys.finalized()
return fsys, potential
def make_zigzag_lead(W, pot=0.0, init = (0,0), a0 = a0, t=2.8, direction = -1):
glat = kwant.lattice.honeycomb(a0, norbs=1)
sym = kwant.TranslationalSymmetry((direction*a0, 0))
sym.add_site_family(glat.sublattices[0], other_vectors=[(-1, 2)])
sym.add_site_family(glat.sublattices[1], other_vectors=[(-1, 2)])
lead = kwant.Builder(sym)
def lead_shape(pos):
x, y = pos
return init[1] <= y <= W+init[1]
def hopping(sitei, sitej, p):
# FROM sitej TO sitei
xi, yi = sitei.pos
xj, yj = sitej.pos
return -t * exp(1j * prefac * p.B * (xi - xj) * (yi + yj) * 0.5*1e-18)
lead[glat.shape(lead_shape, init)] = pot
lead[glat.neighbors()] = hopping
lead.eradicate_dangling()
return lead
leng = 60
wid = 80
B = 0
w = 5.0
pottype = "linear"
V0 = 0.4
E = 0
xcol_realpos = leng - 5 # x where I want the Husimi
par = types.SimpleNamespace(B = B)
pnsys, pot = make_graphene_pnjunction(leng, wid, V0=V0, w=w, pottype=pottype)
pnwf = kwant.solvers.default.wave_function(pnsys, energy = E, args=[par])(0)
pnJ = kwant.operator.Current(pnsys)
M = pnwf.shape[0]
A = pnsys.sites[0].family
def family_shape(i):
site = pnsys.sites[i]
return ('p', 3, 180) if site.family == A else ('p', 3, 0)
# Finally plot stuff:
plt.figure()
for (k, m) in enumerate([0,8]):
wf = np.abs(pnwf[m, :])**2
ax = plt.subplot(1,2,k+1)
kwant.plot(pnsys, site_color=wf, site_symbol=family_shape,
site_size=0.5, hop_lw=0, cmap='BuPu_r', ax = ax)
plt.suptitle("using plot")
plt.figure()
for (k, m) in enumerate([0,8]):
wf = np.abs(pnwf[m, :])**2
ax = plt.subplot(1,2,k+1)
kwant.plotter.map(pnsys, wf, cmap="BuPu_r", ax=ax, oversampling=3)
plt.suptitle("using map")
plt.figure()
for (k, m) in enumerate([0,8]):
wf = pnwf[m, :]
ax = plt.subplot(1,2,k+1)
kwant.plotter.current(pnsys, pnJ(wf, args=[par]), cmap="BuPu_r", ax=ax,
linecolor = (0,0,0, 0.25), min_linewidth=0.5, max_linewidth=1.5)
plt.suptitle("using current")