Dear Pablo,
Thank You very much for Your really helpful replies!
I looked into Your advices and solution, I understood the code, and now some
new issues appear.
Thanks to You, I can see now a beautiful QHE in the energy gap, but I am worry
about the large noises in transverse conductivity beyond the energy gap (I
suppose in this region it should be zero). The second thing is that the energy
gap is not in the zero-energy, but around 4t. Thus, I want to check how the
band structure looks like for this system with given parametres.
According to the KWANT documentation, I found a method to calculate the band
structure: kwant.physics.Bands() but it works only for infinite system.
Therefore, I've created two systems:
1. with leads - to calculate a band structure
2. closed (the previous one) - to calculate the components of conductvity
tensor with local vectors defined away from the edges (as You adviced).
But here I encounter a problem with the first one: despite I know that the
system do has attached leads (so it is infinite), I get an error that the
method kwant.plotter.bands() is 'Expecting an instance of InfiniteSystem.'
The second question is (assuming that this solution would works): are the
results obtained for this two systems can be related with each other? Means,
can I assume that the band structure plotted for the first system will
corresponds to the second system, where I obtained the result for conductivites?
I think also about a different approach: to create a system with four leads and
to calculate (longitudinal and transverse) conductivities as a transmissions
throught a proper leads using the scattering matrix. Maybe that solution could
be more convenient for this problem.
Best regards,
Anna
Below I attach the whole code.
#------------------------------------------------------------------
import kwant
from matplotlib import pyplot
import tinyarray
import numpy as np
# we define the identity matrix and Pauli matrices
sigma_0 = tinyarray.array([[1,0], [0,1]])
sigma_x = tinyarray.array([[0,1], [1,0]])
sigma_y = tinyarray.array([[0,-1j], [1j,0]])
sigma_z = tinyarray.array([[1,0], [0,-1]])
def make_system(t=1, alpha=1, e_z=0.5, W=30, L=30, a=1): # this parameters can
be overwritten in the 'main()' section
lat = kwant.lattice.square(norbs=2) # norbs - the number of orbitals per
site
syst = kwant.Builder()
# Zeeman energy adds to the on-site term
syst[(lat(x,y) for x in range(L) for y in range(W))] = 4*t*sigma_0 +
e_z*sigma_z
# hoppings in x-direction; (a,0) means hopping form (i,j) to (i+1,j)
syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_z+
1j*alpha*sigma_y/(2*a)
# hoppings in y-direction; (0,a) means hopping form (i,j) to (i,j+1)
syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_z +
1j*alpha*sigma_x/(2*a)
return syst
def make_system_leads(t=1, alpha=1, e_z=0.5, W=30, L=30, a=1): # this
parameters can be overwritten in the 'main()' section
lat = kwant.lattice.square(norbs=2) # norbs - the number of orbitals per
site
syst = kwant.Builder()
# Zeeman energy adds to the on-site term
syst[(lat(x,y) for x in range(L) for y in range(W))] = 4*t*sigma_0 +
e_z*sigma_z
# hoppings in x-direction; (a,0) means hopping form (i,j) to (i+1,j)
syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_z+
1j*alpha*sigma_y/(2*a)
# hoppings in y-direction; (0,a) means hopping form (i,j) to (i,j+1)
syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_z +
1j*alpha*sigma_x/(2*a)
# create the leads (left and right)
lead = kwant.Builder(kwant.TranslationalSymmetry((-a,0)))
lead[(lat(0,j) for j in range(W))] = 4*t*sigma_0 + e_z*sigma_z
lead[kwant.builder.HoppingKind((1,0), lat, lat)] = -t*sigma_0 -
1j*alpha*sigma_y/(2*a) # hoppings in x-direction
lead[kwant.builder.HoppingKind((0,1), lat, lat)] = -t*sigma_0 +
1j*alpha*sigma_x/(2*a) # hoppings in y-direction
# attach the leads (left and right)
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
# create the leads (top and bottom)
#lead2 = kwant.Builder(kwant.TranslationalSymmetry((0,-a)))
#lead2[(lat(i,0) for i in range(L))] = 4*t*sigma_0 + e_z*sigma_z
#lead2[kwant.builder.HoppingKind((1,0), lat, lat)] = -t*sigma_0 -
1j*alpha*sigma_y/(2*a) # hoppings in x-direction
#lead2[kwant.builder.HoppingKind((0,1), lat, lat)] = -t*sigma_0 +
1j*alpha*sigma_x/(2*a) # hoppings in y-direction
# attach the leads (top and bottom)
#syst.attach_lead(lead2)
#syst.attach_lead(lead2.reversed())
return syst
def plot_dos_and_curves(dos, labels_to_data):
pyplot.figure()
pyplot.fill_between(dos[0], dos[1], label="DoS [a.u.]",
alpha=0.5, color='gray')
for label, (x, y) in labels_to_data:
pyplot.plot(x, y, label=label, linewidth=2)
pyplot.legend(loc=2, framealpha=0.5)
pyplot.ylim(-1, 3) #!
pyplot.xlabel("energy [t]")
pyplot.ylabel("$\sigma [e^2/h]$")
pyplot.show()
def main(random_vecs=True):
num_moments = 400 # used in the method kwant.kpm.SpectralDensity
# 'num_moments' is the number of moments, order of the KPM expansion.
Mutually exclusive with 'energy_resolution'
# default num_moments = 100
# dimensions and parameters for the both systems
W = 30
L = 30
alpha = 1
e_z = 0.5
# calculate the band structure (infinite system: lead on the left and on
the right hand side)
syst_infty = make_system_leads(W=W, L=L, alpha=alpha, e_z=e_z)
kwant.plot(syst_infty); # check that the system looks as intended (with
the left and right leads attached)
syst_infty = syst_infty.finalized()
kwant.plotter.bands(syst_infty) # HERE I GET AN ERROR: 'Expecting an
instance of InfiniteSystem.'
# the closed system to calculate conductivities using the method:
kwant.kpm.conductivity()
syst = make_system(W=W, L=L, alpha=alpha, e_z=e_z)
syst = syst.finalized()
fam = syst.sites[0].family # assuming all sites from the same
(sub)lattice
area_per_orb = np.cross(*fam.prim_vecs) / fam.norbs # area that each of
our vectors covers
#----------------------------------------------------------
if random_vecs:
edge_width = min(L, W) / 4
def center(s):
return (edge_width <= s.pos[0] < L - edge_width) and (edge_width <=
s.pos[1] < W - edge_width)
center_vectors = kwant.kpm.RandomVectors(syst, where=center)
# here kwant.kpm.LocalVectors can also be used, but will produce much
more vectors
one_vector = next(center_vectors)
num_vectors = 10 # len(center_vectors)
else: # use local vectors
edge_width = min(L, W) / 4
def center(s):
return (edge_width <= s.pos[0] < L - edge_width) and (edge_width <=
s.pos[1] < W - edge_width)
center_vectors = np.array(list(kwant.kpm.LocalVectors(syst,
where=center)))
one_vector = center_vectors[0]
num_vectors = len(center_vectors) # len(x) - length of the object 'x'
#----------------------------------------------------------
# we need to normalize the results by the area that each of our vectors
covers
norm = area_per_orb * np.linalg.norm(one_vector) ** 2
# check the area that the vectors cover
kwant.plot(syst, site_color=np.abs(one_vector[::2]) +
np.abs(one_vector[1::2]), cmap='Reds');
# object that represents the density of states for the system
spectrum = kwant.kpm.SpectralDensity(syst, rng=0, num_moments=num_moments,
num_vectors=num_vectors,
vector_factory=center_vectors)
# vector_factory - optional; it should contain (or yield) vectors of the
size of the system
# (here we use only the vector in the center of the system, thus we need to
define this argument;
# otherwise the KWNAT will randomize vectors for whole sample area)
energies, densities = spectrum()
# component 'xx'
cond_xx = kwant.kpm.conductivity(syst, alpha='x', beta='x',
num_moments=num_moments,
rng=0, num_vectors=num_vectors,
vector_factory=center_vectors)
# component 'xy'
cond_xy = kwant.kpm.conductivity(syst, alpha='x', beta='y',
num_moments=num_moments,
rng=0,
num_vectors=num_vectors,vector_factory=center_vectors)
energies = cond_xx.energies
cond_array_xx = np.array([cond_xx(e, temperature=0.01) for e in energies])
/ norm
cond_array_xy = np.array([cond_xy(e, temperature=0.01) for e in energies])
/ norm
scale_to_plot = 1 / (edge_width)
plot_dos_and_curves(
(spectrum.energies, spectrum.densities / norm),
[(r'Longitudinal conductivity $\sigma_{xx}$ (scaled)',
(energies, cond_array_xx.real * scale_to_plot)),
(r'Hall conductivity $\sigma_{xy}$',
(energies, cond_array_xy.real))],
)
# standard Python construct to execute 'main' if the program is used as a script
if __name__ == '__main__':
main()
