plotting band structure of bulk 2D material along high symmetry lines in the Brillouin zone
Dear kwant developers and users, I want to calculate and plot the band structure of bulk 2D material along some highsymmetry lines in the Brillouin zone (e.g. Gamma>K>MGamma), 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://kwantdiscuss.kwantproject.narkive.com/2fT27IJh/kwantforbulksyst.... 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) > 1e7): 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 highsymmetry lines in the Brillouin zone (e.g. Gamma>K>MGamma), 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://kwantdiscuss.kwantproject.narkive.com/2fT27IJh/kwantforbulksyst....
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
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) > 1e7): 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 highsymmetry lines in the Brillouin zone (e.g. Gamma>K>MGamma), 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://kwantdiscuss.kwantproject.narkive.com/2fT27IJh/kwantforbulksyst....
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
participants (2)

Abbout Adel

kuangyia lee