Hello Kwant community, I have played with a Weyl semimetal system for a while, using Kwant's conductance features (with leads). And I do get the quantized Hall response under magnetic field as reported in this work [PRL 125, 036602]. Now with exactly the same system (size, Fermi energy, hopping, etc.), I remove the leads and try the kpm.conductivity features. However, there is no sign of quantization at all, as far as I've tried with LocalVectors or RandomVectors and tuning num_moments or so. To make sure of my use of kpm.conductivity, I tried another known unquantized Hall effect (sigma_xy) in this system without magnetic field and it looked not unreasonable; I also changed the system to a simple 2D QHE one and the same code worked well to see quantization. This seems rather unclear to me. Could you please help with how should I correctly calculate the conductivity in Kwant here? The code below plots sigma_xz and sigma_zx with respect to a few inverse magnetic fields. Thanks! from cmath import exp from math import cos, sin, pi import numpy as np import time import sys import kwant # For plotting from matplotlib import pyplot as plt # For matrix support import tinyarray # define Pauli-matrices for convenience 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]]) Lx=25 Ly=20 Lz=25 D=0.02 D1=D D2=4*D A=1 M=0.3 kw=pi/2.0 E_F=2*D2*(1-cos(kw)) # Fermi energy norbs=2 def make_system(a=1, t=1): lat = kwant.lattice.cubic(a, norbs=norbs) syst = kwant.Builder() def fluxx(site1, site2, B): pos = [(site1.pos[i]+site2.pos[i])/2.0 for i in range(3)] return exp(-1j * (B * pos[2])) def hopx(site1, site2, B): return fluxx(site1, site2, B) * (-D2 * sigma_0 + 1j * A/2 * sigma_x + M * sigma_z) def hopy(site1, site2): return (-D1 * sigma_0 + 1j * A/2 * sigma_y + M * sigma_z) def hopz(site1, site2): return (-D2 * sigma_0 + M * sigma_z) def onsite(site): return (2 * D1 + 2 * D2 * 2) * sigma_0 + 2 * M * ((1 - cos(kw)) - 3) * sigma_z syst[(lat(x, y, z) for x in range(Lx) for y in range(Ly) for z in range(Lz))] = onsite # hoppings in x-direction syst[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hopx # hoppings in y-directions syst[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hopy # hoppings in z-directions syst[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hopz return lat, syst.finalized() 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$') plt.show() plt.clf() def cond(lat, syst, random_vecs_flag=False): center_width = min(Lx,Ly,Lz)//20 center_pos = [Lx//2,Ly//2,Lz//2] def center(s): return (np.abs(s.pos[0]-center_pos[0]) < center_width) and (np.abs(s.pos[1]-center_pos[1]) < center_width) and (np.abs(s.pos[2]-center_pos[2]) < center_width) num_mmts = 200; if random_vecs_flag: #use RandomVecs center_vecs = kwant.kpm.RandomVectors(syst, where=center) one_vec = next(center_vecs) num_vecs = 20 else: # use LocalVecs center_vecs = np.array(list(kwant.kpm.LocalVectors(syst, where=center))) one_vec = center_vecs[0] num_vecs = len(center_vecs) area_per_orb = np.abs(np.linalg.det(lat.prim_vecs)) norm = area_per_orb * np.linalg.norm(one_vec) ** 2 / norbs Bfields = [] Binvs = np.arange(2.4, 65.2, 30.5) for BInv in Binvs: Bfields.append(1/BInv) params = dict(B=0) cond_array_zx = [] cond_array_xz = [] for B in Bfields: params['B'] = B cond_zx = kwant.kpm.conductivity(syst, params=params, alpha='z', beta='x', mean=True, rng=0, num_moments=num_mmts, num_vectors=num_vecs, vector_factory=center_vecs) cond_xz = kwant.kpm.conductivity(syst, params=params, alpha='x', beta='z', mean=True, rng=0, num_moments=num_mmts, num_vectors=num_vecs, vector_factory=center_vecs) cond_array_zx.append(cond_zx(mu=E_F, temperature=0.0)/ norm) cond_array_xz.append(cond_xz(mu=E_F, temperature=0.0)/ norm) cond_array_zx = np.array(cond_array_zx) * Ly # 3D conductivity * Ly gives the sheet 2D-like conductivity supposed to be quantized cond_array_xz = np.array(cond_array_xz) * Ly print(cond_array_zx) print(cond_array_xz) print('Calculating duration (s) until finalising system: ', time.time() - global_start_time) plot_curves([(r'$\sigma_{zx}$', (Binvs, cond_array_zx)), (r'$\sigma_{xz}$', (Binvs, cond_array_xz))]) def main(): lat, syst = make_system() # Check that the system looks as intended. #kwant.plot(syst, file="lattice_fig.pdf") cond(lat, syst) if __name__ == '__main__': global_start_time = time.time() main()