Hello,
I encountered an issue when dealing with a semi-infinite ribbon: when defining a 1D translationnal invariance, I end up with a unit cell with unwanted and unexpected additionnal degrees of freedom.
The unit cell have additionnal segments of atoms (depending on the translationnal invariance vector). This additionnal segment is not visible with a "kwant.plot()" but appears in the Hilbert space dimensionnality.
It then create unexpected states close to zero energy which are not visible on the band structure but are visible in the DOS.
I copy the script at the end of the message, if you want to have a look. This is a Haldane model ribbon with armchair edges.
I tried to change the translationnal invariance vector but I didn't manage to get around it.
Do you have a solution ? Should I use wraparound ?
My end goal is to have a large enough system with a kind of periodic boundary condition in one direction, but still preserving a plotable dispersion relation, so I though that using 1D translationnal invariance with a big unit cell was a good idea.
Thank you in advance,
Alexandre BERNARD
---------------------------------------------------- the script ---------------------------------------------------------------
from types import SimpleNamespace
import matplotlib
import math
from matplotlib import pyplot
from mpl_toolkits import mplot3d
import numpy as np
import scipy.linalg as la
import random
import kwant
def random_uni(param):
return param*(random.random() - 0.5)*2 # returns a random number between +param and -param
def onsite(site, p):
# Here we affect a potential whose sign depends on which sublattice the atom lives in (breaking sublattice symmetry)
if site.family == A: # sublattice A
return p.m + random_uni(p.rand_os)
else: # sublattice B
return -p.m + random_uni(p.rand_os)
def nn_hopping(site1, site2, params):
return params.t_1
def nnn_hopping(site1, site2, params):
return (math.cos(params.phi) + 1j*math.sin(params.phi)) * (params.t_2)
def ribbon_shape_armchair(pos):
return (-1 <= pos[0] < w)
def hopping_color(site1, site2):
if site1.family == A and site2.family == A:
return 'b'
elif site1.family == B and site2.family == B:
return 'r'
else:
return '0.5'
def site_color(site):
if site.family == A:
return 'b'
else:
return 'r'
np.random.seed(271828) # Fixing the RNG seed for reproducibility
# Graphene lattice "pointing up", with two sublattices
graphene = kwant.lattice.general([[1, 0], [1/2, np.sqrt(3)/2]], # Here are the vectors of the Bravais lattice
[[0, 0], [0, 1/np.sqrt(3)]]) # Here are the positions of the atoms in the unit cell
A, B = graphene.sublattices
w = 5 # width of the ribbon
p = SimpleNamespace(t_1=1.0, t_2=0.1, m=0.0, phi=np.pi/2, rand_os=0.)
# Creation of the list of (oriented) bonds that will have the complex hopping term
nnn_hoppings_A = (((-1, 0), A, A), ((0, 1), A, A), ((1, -1), A, A)) # Bonds within the A sublattice
nnn_hoppings_B = (((1, 0), B, B), ((0, -1), B, B), ((-1, 1), B, B)) # Bonds within the B sublattice
nnn_hoppings = nnn_hoppings_A + nnn_hoppings_B
haldane_armchair = kwant.Builder(kwant.TranslationalSymmetry((0, 1*np.sqrt(3))))
haldane_armchair[graphene.shape(ribbon_shape_armchair, (0, 0))] = onsite
haldane_armchair[graphene.neighbors(1)] = nn_hopping
haldane_armchair[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings]] = nnn_hopping
kwant.plot(haldane_armchair, site_color=site_color, hop_color=hopping_color, fig_size=(6, 9));
haldane_f2 = haldane_armchair.finalized()
ham = haldane_f2.hamiltonian_submatrix(args = [p])
eigenVal, eigenVec = la.eigh(ham)
nb_eVal = len(eigenVal[:])
print("Number of eigenvalues (i.e. unit cell DOF): ", nb_eVal)
pyplot.hist(eigenVal[:], 300)
pyplot.show()
pyplot.hist(eigenVal[:], 100, cumulative=True)
pyplot.show()
kwant.plotter.bands(haldane_armchair.finalized(), args=[SimpleNamespace(t_1=p.t_1, t_2=p.t_2, m=p.m, phi=p.phi, rand_os=p.rand_os)], fig_size=(12, 8));
E_F = 0.0
for j in range(nb_eVal):
if eigenVal[j] > E_F:
break
print("eigenVal[", j, "] = ", eigenVal[j])
pyplot.figure(figsize=(w,w/2))
for i, s in enumerate(haldane_f2.sites):
x,y=s.pos
prob = abs(eigenVec[i,j])**2
pyplot.scatter(x,y,c=prob,vmin=0,vmax=5e-3,s=60)
pyplot.show()
