Dear kwant developers and users,
I want to calculate and plot the band structure of bulk 2D material along some high-symmetry lines in the Brillouin zone (e.g. Gamma-->K-->M--Gamma), i.e. the 1D band structure.
I found that someone already posted a relevant issue and some useful solutions were also given. For example, see https://kwant-discuss.kwant-project.narkive.com/2fT27IJh/kwant-for-bulk-syst....
However, these solutions are only for calculating and plotting the full 2D band structure in the Brillouin zone.
Can anyone give some useful solutions to my concern? Your help is much appreciated.
Thank you and best regards, Kuangyia Lee
Dear Lee,
All what you need is to express the equation of your path: I mean ky as a function of kx on this path. Please, find below an example about Graphene band structure. I compare it to the analytical result found in literature.
Just be careful, when the lattice vectors are not orthogonal. You need to do some transformation to get the vectors in the reciprocal lattice.
I wrote this a while ago, so please check it if there are no mistakes in it. It works with kwant 1.3
I hope this helps Adel ###############################
import numpy as np from numpy import pi,sqrt,arccos import kwant from kwant.wraparound import wraparound,plot_2d_bands from matplotlib import pyplot import scipy
def make_graphene(): """Create a builder for a periodic infinite sheet of graphene.""" lat = kwant.lattice.honeycomb() sym = kwant.TranslationalSymmetry(lat.vec((1, 0)), lat.vec((0, 1))) sys = kwant.Builder(sym) sys[lat.shape(lambda p: True, (0, 0))] = 0 sys[lat.neighbors()] = 1 return sys
sys=make_graphene() syst = wraparound(sys).finalized() kwant.plot(sys) plot_2d_bands(syst)
# calculate the reciprocal symmetry vectors: these # are the columns of 'A' sym=sys.symmetry B = np.array(sym.periods).T A = B.dot(np.linalg.inv(B.T.dot(B))) # transforms a momentum in the basis (kx, ky) to the basis of # reciprocal symmetry vectors def momentum_to_lattice(k): k, residuals = scipy.linalg.lstsq(A, k)[:2] if np.any(abs(residuals) > 1e-7): raise RuntimeError("Requested momentum doesn't correspond" " to any lattice momentum.") return k
def H_k(kx, ky): k = (kx, ky) k_prime = momentum_to_lattice(k) return syst.hamiltonian_submatrix(args=k_prime)
# This function is the exact result taken from literature def function(kx,ky): return sqrt(3+2*np.cos(kx)+4*np.cos(kx/2.) *np.cos(sqrt(3)/2*ky))
# Here I am comparing the numerical result with the analytical result. with kx varying from -pi to pi k_val=np.linspace(-pi,pi,61) ky=2*pi/sqrt(3) resultat1=[abs(H_k(kx,ky)[0]) for kx in k_val] resultat2=[abs(function(kx,ky)) for kx in k_val]
pyplot.plot(k_val,resultat1,'r-') pyplot.plot(k_val,resultat2,'go') pyplot.show()
# # # $ # E=\pm t \sqrt{3+2 \cos(\sqrt{3}k_x a)+a \cos(\frac{\sqrt{3}}{2}k_x a) \cos(\frac{3}{2} k_y a)} # $ # # with $a=1/\sqrt{3}$ # # The point $K$ is: $K(\frac{4\pi}{3\sqrt{3}a},0) $ # The point $M$ is: $M(\frac{\pi}{\sqrt{3} a},\frac{\pi}{3 a})$
# In[ ]:
def eigenvalues(kx,ky): return np.sort(np.linalg.eigvalsh(H_k(kx, ky)))
nb=30 d1=4*pi/3 # this is the distance when you go from Gamma to K d2=2*pi/3 # this is the distance when you go from K to M
Kx1= np.linspace(0,4*pi/3,nb) Result1=np.array([eigenvalues(kx,0) for kx in Kx1])
Kx2=[d1 + sqrt((kx - (4*pi/3))**2 + 3*(kx - 4*pi/3)**2) for kx in np.linspace(pi,4*pi/3,nb)] #distance Gamma K is added here Result2=np.array([eigenvalues(kx,-sqrt(3)*(kx - 4.* pi/3)) for kx in np.linspace(pi,4*pi/3,nb)])
Kx3=[d1 + d2 + sqrt((x - pi)**2 + (1/sqrt(3)* x - pi/sqrt(3))**2) for x in np.linspace(0,pi,nb)]##3distance Gamma K M is addded here Result3=np.array([eigenvalues(kx, 1./sqrt(3)* kx) for kx in np.linspace(0,pi,nb)])
pyplot.scatter(Kx1,Result1[:,0],color='r') pyplot.scatter(Kx2,Result2[:,0],color='g') pyplot.scatter(Kx3,Result3[:,0],color='c')
pyplot.scatter(Kx1,Result1[:,1],color='r') pyplot.scatter(Kx2,Result2[:,1],color='g') pyplot.scatter(Kx3,Result3[:,1],color='c')
pyplot.show()
On Wed, Oct 24, 2018 at 7:53 PM kuangyia lee kuangyia.lee@gmail.com wrote:
Dear kwant developers and users,
I want to calculate and plot the band structure of bulk 2D material along some high-symmetry lines in the Brillouin zone (e.g. Gamma-->K-->M--Gamma), i.e. the 1D band structure.
I found that someone already posted a relevant issue and some useful solutions were also given. For example, see https://kwant-discuss.kwant-project.narkive.com/2fT27IJh/kwant-for-bulk-syst....
However, these solutions are only for calculating and plotting the full 2D band structure in the Brillouin zone.
Can anyone give some useful solutions to my concern? Your help is much appreciated.
Thank you and best regards, Kuangyia Lee
Dear Adel,
Thank you for the tutorial script.
Best, Kuangyia
On Wed, Oct 24, 2018 at 11:33 PM Abbout Adel abbout.adel@gmail.com wrote:
Dear Lee,
All what you need is to express the equation of your path: I mean ky as a function of kx on this path. Please, find below an example about Graphene band structure. I compare it to the analytical result found in literature.
Just be careful, when the lattice vectors are not orthogonal. You need to do some transformation to get the vectors in the reciprocal lattice.
I wrote this a while ago, so please check it if there are no mistakes in it. It works with kwant 1.3
I hope this helps Adel ###############################
import numpy as np from numpy import pi,sqrt,arccos import kwant from kwant.wraparound import wraparound,plot_2d_bands from matplotlib import pyplot import scipy
def make_graphene(): """Create a builder for a periodic infinite sheet of graphene.""" lat = kwant.lattice.honeycomb() sym = kwant.TranslationalSymmetry(lat.vec((1, 0)), lat.vec((0, 1))) sys = kwant.Builder(sym) sys[lat.shape(lambda p: True, (0, 0))] = 0 sys[lat.neighbors()] = 1 return sys
sys=make_graphene() syst = wraparound(sys).finalized() kwant.plot(sys) plot_2d_bands(syst)
# calculate the reciprocal symmetry vectors: these # are the columns of 'A' sym=sys.symmetry B = np.array(sym.periods).T A = B.dot(np.linalg.inv(B.T.dot(B))) # transforms a momentum in the basis (kx, ky) to the basis of # reciprocal symmetry vectors def momentum_to_lattice(k): k, residuals = scipy.linalg.lstsq(A, k)[:2] if np.any(abs(residuals) > 1e-7): raise RuntimeError("Requested momentum doesn't correspond" " to any lattice momentum.") return k
def H_k(kx, ky): k = (kx, ky) k_prime = momentum_to_lattice(k) return syst.hamiltonian_submatrix(args=k_prime)
# This function is the exact result taken from literature def function(kx,ky): return sqrt(3+2*np.cos(kx)+4*np.cos(kx/2.) *np.cos(sqrt(3)/2*ky))
# Here I am comparing the numerical result with the analytical result. with kx varying from -pi to pi k_val=np.linspace(-pi,pi,61) ky=2*pi/sqrt(3) resultat1=[abs(H_k(kx,ky)[0]) for kx in k_val] resultat2=[abs(function(kx,ky)) for kx in k_val]
pyplot.plot(k_val,resultat1,'r-') pyplot.plot(k_val,resultat2,'go') pyplot.show()
# # # $ # E=\pm t \sqrt{3+2 \cos(\sqrt{3}k_x a)+a \cos(\frac{\sqrt{3}}{2}k_x a) \cos(\frac{3}{2} k_y a)} # $ # # with $a=1/\sqrt{3}$ # # The point $K$ is: $K(\frac{4\pi}{3\sqrt{3}a},0) $ # The point $M$ is: $M(\frac{\pi}{\sqrt{3} a},\frac{\pi}{3 a})$
# In[ ]:
def eigenvalues(kx,ky): return np.sort(np.linalg.eigvalsh(H_k(kx, ky)))
nb=30 d1=4*pi/3 # this is the distance when you go from Gamma to K d2=2*pi/3 # this is the distance when you go from K to M
Kx1= np.linspace(0,4*pi/3,nb) Result1=np.array([eigenvalues(kx,0) for kx in Kx1])
Kx2=[d1 + sqrt((kx - (4*pi/3))**2 + 3*(kx - 4*pi/3)**2) for kx in np.linspace(pi,4*pi/3,nb)] #distance Gamma K is added here Result2=np.array([eigenvalues(kx,-sqrt(3)*(kx - 4.* pi/3)) for kx in np.linspace(pi,4*pi/3,nb)])
Kx3=[d1 + d2 + sqrt((x - pi)**2 + (1/sqrt(3)* x - pi/sqrt(3))**2) for x in np.linspace(0,pi,nb)]##3distance Gamma K M is addded here Result3=np.array([eigenvalues(kx, 1./sqrt(3)* kx) for kx in np.linspace(0,pi,nb)])
pyplot.scatter(Kx1,Result1[:,0],color='r') pyplot.scatter(Kx2,Result2[:,0],color='g') pyplot.scatter(Kx3,Result3[:,0],color='c')
pyplot.scatter(Kx1,Result1[:,1],color='r') pyplot.scatter(Kx2,Result2[:,1],color='g') pyplot.scatter(Kx3,Result3[:,1],color='c')
pyplot.show()
On Wed, Oct 24, 2018 at 7:53 PM kuangyia lee kuangyia.lee@gmail.com wrote:
Dear kwant developers and users,
I want to calculate and plot the band structure of bulk 2D material along some high-symmetry lines in the Brillouin zone (e.g. Gamma-->K-->M--Gamma), i.e. the 1D band structure.
I found that someone already posted a relevant issue and some useful solutions were also given. For example, see https://kwant-discuss.kwant-project.narkive.com/2fT27IJh/kwant-for-bulk-syst....
However, these solutions are only for calculating and plotting the full 2D band structure in the Brillouin zone.
Can anyone give some useful solutions to my concern? Your help is much appreciated.
Thank you and best regards, Kuangyia Lee
-- Abbout Adel