Hello Kwant community, I am using the KPM method to calculate conductivities in the simplest QHE system -- square lattice nearest hopping gas. The goal is to see well-quantized plateaus up to n=5. I use local vectors at the center of the system. But it seems not available at all as far as I've tried increasing system size (like 80*80) and number of moments (like 500). I doubt if it really needs so demanding computational resources in order to get n=5. Hopefully I missed something basic here. Thank you! from cmath import exp import numpy as np import time import sys import kwant #from matplotlib import pyplot from matplotlib import pyplot as plt L_sys = [50, 50] E_F=0.4 def make_system(a=1, t=1): lat = kwant.lattice.square(a, norbs=1) syst = kwant.Builder() def fluxx(site1, site2, Bvec): pos = [(site1.pos[i]+site2.pos[i]-L_sys[i])/2.0 for i in range(2)] return exp(-1j * Bvec * pos[1] ) def fluxy(site1, site2, Bvec): pos = [(site1.pos[i]+site2.pos[i]-L_sys[i])/2.0 for i in range(2)] return 1 #exp(-1j * (-Bvec * pos[0])) def hopx(site1, site2, Bvec): # The magnetic field is controlled by the parameter Bvec return fluxx(site1, site2, Bvec) * t def hopy(site1, site2, Bvec): # The magnetic field is controlled by the parameter Bvec return fluxy(site1, site2, Bvec) * t def onsite(site): return 4*t #### Define the scattering region. #### syst[(lat(x, y) for x in range(L_sys[0]) for y in range(L_sys[1]))] = onsite # hoppings in x-direction syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx # hoppings in y-directions syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy # Finalize the system. syst = syst.finalized() return lat, syst def cond(lat, syst, random_vecs_flag=False): center_pos = [x/2 for x in L_sys] center0 = lambda s: np.linalg.norm(s.pos-center_pos) < 2 num_mmts = 800; if random_vecs_flag: center_vecs = kwant.kpm.RandomVectors(syst, where=center0) one_vec = next(center_vecs) num_vecs = 10 # len(center_vecs) else: center_vecs = np.array(list(kwant.kpm.LocalVectors(syst, where=center0))) one_vec = center_vecs[0] num_vecs = len(center_vecs) area_per_orb = np.abs(np.cross(*lat.prim_vecs)) norm = area_per_orb * np.linalg.norm(one_vec) ** 2 Binvfields = np.arange(3.01, 26.2, 1) params = dict(Bvec=0) dataxx = [] dataxy = [] for Binv in Binvfields: params['Bvec'] = 1/Binv cond_xx = kwant.kpm.conductivity(syst, params=params, alpha='x', beta='x', mean=True, num_moments=num_mmts, num_vectors=num_vecs, vector_factory=center_vecs) cond_xy = kwant.kpm.conductivity(syst, params=params, alpha='x', beta='y', mean=True, num_moments=num_mmts, num_vectors=num_vecs, vector_factory=center_vecs) dataxx.append(5*cond_xx(mu=E_F, temperature=0.0)/norm) dataxy.append(-cond_xy(mu=E_F, temperature=0.0)/norm) print(dataxx) print(dataxy) plot_curves( [ (r'$\sigma_{xx}\times5$',(Binvfields,dataxx)), (r'$\sigma_{xy}$',(Binvfields,dataxy)), ] ) def plot_curves(labels_to_data): plt.figure(figsize=(12,10)) for label, (x, y) in labels_to_data: plt.plot(x, y, label=label, linewidth=2) plt.legend(loc=2, framealpha=0.5) plt.xlabel("1/B") plt.ylabel(r'$\sigma [e^2/h]$') plt.show() plt.clf() def main(): lat, syst = make_system() cond(lat, syst) if __name__ == '__main__': main()