Dear Anna,
I modified your Hamiltonian and exaggerated the value of the parameters in
order to have a larger gap.
If you want to use the values that you had the first time, your gap will be
much smaller. In that case, you need more energy resolution in the KPM
expansion, and that will slow down things a lot (just beware). Likely you will
need a larger sample when the gap is small.
```
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.8, alpha=0.1, e_z=0.0018, W=30, L=30, a=1):
lat = kwant.lattice.square(norbs=2) # norbs - the number of orbitals per
site
syst = kwant.Builder()
syst[(lat(x,y) for x in range(L) for y in range(W))] = 4*t*sigma_0 +
e_z*sigma_z # Zeeman energy adds to the on-site term
syst[kwant.builder.HoppingKind((a,0), lat, lat)] = -t*sigma_z+
1j*alpha*sigma_y/(2*a) #
# hoppings in x-direction; (a,0) means hopping form (i,j) to (i+1,j)
syst[kwant.builder.HoppingKind((0,a), lat, lat)] = -t*sigma_z +
1j*alpha*sigma_x/(2*a) #
# hoppings in y-direction; (0,a) means hopping form (i,j) to (i,j+1)
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
W=30
L=30
syst = make_system(W=W, L=L, alpha=1, e_z=0.5)
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
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)
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)
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');
spectrum = kwant.kpm.SpectralDensity(syst, rng=0,num_moments=num_moments,
num_vectors=num_vectors, vector_factory=center_vectors) # object
# that represents the density of states for the system
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))],
)
```
