Dear developer, I was trying to calculate the conductivity of honeycomb lattice, where I have included the Haldane term along with the spin-orbit Rashba term (Kane-Mele model). I was following the tutorial provided in chapter 2.9 (https://kwant-project.org/doc/1/tutorial/kpm). The code shows an error when I include the spin in the system. My code is attached below, from matplotlib import pyplot as plt import tinyarray as tn import kwant import numpy as np t = 1 t2 = 0.3 sgm_0=tn.array([[1,0],[0,1]]) sgm_x=tn.array([[0,1],[1,0]]) sgm_y=tn.array([[0,-1j],[1j,0]]) sgm_z=tn.array([[1,0],[0,-1]]) def make_syst_topo(l=20, a=1, m=0): alpha = 0.01 syst = kwant.Builder() lat = kwant.lattice.honeycomb(a, norbs=1, name=['a', 'b']) a = lat.a b = lat.b def rectangle(pos): x, y = pos return - l<= x <= +l and - l<= y <= +l def rashba(site_i,site_j): d_ij = site_i.pos - site_j.pos return -1j*alpha*(sgm_x*d_ij[1]-sgm_y*d_ij[0]) def nn_hop(site_i, site_j): return rashba(site_i, site_j) - t * sgm_0 syst[lat.shape(rectangle, (0, 0))] = 0.0*sgm_0 syst[kwant.builder.HoppingKind((0, 0), a, b)] = nn_hop syst[kwant.builder.HoppingKind((0, 1), a, b)] = nn_hop syst[kwant.builder.HoppingKind((-1, 1), a, b)] = nn_hop syst[lat.a.neighbors()] = 1j * t2 * sgm_z syst[lat.b.neighbors()] = -1j * t2*sgm_z syst.eradicate_dangling() return lat, syst.finalized() def plot_system(syst): kwant.plot(syst,site_size=0.1, hop_lw=0.01) # Plot fill density of states plus curves on the same axes. def plot_dos_and_curves(dos, labels_to_data): plt.figure(figsize=(10,8)) plt.fill_between(dos[0], dos[1], label="DoS [a.u.]", alpha=0.5, color='gray') for label, (x, y) in labels_to_data: plt.plot(x, y, label=label, linewidth=2) plt.legend(loc='upper center', framealpha=0.5, fontsize=12) plt.xlabel(r'$E/t$'); # plt.ylabel('ylabel') plt.show() plt.clf() def conductivity_example(): # construct the Haldane model lat, fsyst = make_syst_topo() # find 'A' and 'B' sites in the unit cell at the center of the disk where = lambda s: np.linalg.norm(s.pos) < 1 # component 'xx' s_factory = kwant.kpm.LocalVectors(fsyst, where) cond_xx = kwant.kpm.conductivity(fsyst, alpha='x', beta='x', mean=True, num_vectors=None, vector_factory=s_factory) # component 'xx' s_factory = kwant.kpm.LocalVectors(fsyst, where) cond_yy = kwant.kpm.conductivity(fsyst, alpha='y', beta='y', mean=True, num_vectors=None, vector_factory=s_factory) # component 'xy' s_factory = kwant.kpm.LocalVectors(fsyst, where) cond_xy = kwant.kpm.conductivity(fsyst, alpha='x', beta='y', mean=True, num_vectors=None, vector_factory=s_factory) energies_x = cond_xx.energies energies_y = cond_yy.energies cond_array_xx = np.array([cond_xx(e, temperature=0.00) for e in energies_x]) cond_array_xy = np.array([cond_xy(e, temperature=0.00) for e in energies_x]) cond_array_yy = np.array([cond_yy(e, temperature=0.00) for e in energies_y]) # area of the unit cell per site area_per_site = np.abs(np.cross(*lat.prim_vecs)) / len(lat.sublattices) cond_array_xx /= area_per_site cond_array_xy /= area_per_site cond_array_yy /= area_per_site # ldos s_factory = kwant.kpm.LocalVectors(fsyst, where) spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None, vector_factory=s_factory) plot_dos_and_curves( (spectrum.energies, spectrum.densities * 8), [ (r'$\sigma_{xx} / 4$', (energies_x, cond_array_xx / 4)), (r'$\sigma_{xy}$', (energies_x, cond_array_xy)), (r'$\sigma_{yy} / 4$', (energies_y, cond_array_yy / 4))] ) def main(): fsyst = make_syst_topo()[1] plot_system(fsyst) conductivity_example() if __name__ == '__main__': main() I shall look forward to your kind reply.