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))], ) ```