Hi everyone, I would like to reproduce figure 2.c of the manuscript http://doi.org/10.1103/PhysRevB.88.125409 in that they calculate the conductivity of a bilayer zigzag graphene. Even though I get similar transmission values, the steps do not occur at the same values. I'm not sure what my mistake would be. I would appreciate any suggestions. The code is the following: import numpy as np import kwant import matplotlib.pyplot as plt a = 0.142 c0 = 0.5 # to try t = 3.16 t1 = 0.39/3.16 widthh= 23.7 graphene = kwant.lattice.general([(a*np.sqrt(3), 0,0), (a*np.sqrt(3) / 2, 3 / 2 * a,0)], [(0, 0,0), (0, -a,0), (0, a,c0), (0, 0,c0) ]) a1, b1, a2, b2 = graphene.sublattices syst = kwant.Builder() def make_system(length=20, width=widthh): def rectangle(pos): x, y, z = pos return -length / 2 < x < length / 2 and -width / 2 < y < width / 2 and -5 < z < 5 syst[graphene.shape(rectangle, (0, 0,0))] = 0 syst[graphene.neighbors()] = t syst[kwant.builder.HoppingKind((0, 0), b2, a1)] = t1 make_system() symmetry_left = kwant.TranslationalSymmetry((-a*np.sqrt(3), 0, 0)) symmetry_right = kwant.TranslationalSymmetry((a*np.sqrt(3), 0, 0)) lead_left = kwant.Builder(symmetry_left) lead_right = kwant.Builder(symmetry_right) def lead_shape(pos): x, y, z = pos return -widthh / 2 < y < widthh / 2 for lead in [lead_left, lead_right]: lead[graphene.shape(lead_shape, (0, 0, 0))] = 0 lead[graphene.neighbors()] = t lead[kwant.builder.HoppingKind((0, 0), b2, a1)] = t1 syst.attach_lead(lead_left); syst.attach_lead(lead_right); kwant.plot(syst) syst_finalized = syst.finalized() lead_index = 0 lead = syst_finalized.leads[lead_index] # Conductivity calculation data = [] for i in range(51): energias = 0+0.07*3.16*i/50.0+0.000001 smatrix_bilayer = kwant.smatrix(syst_finalized, energy=energias) data.append([energias, smatrix_bilayer.transmission(1, 0)]) data = np.array(data) plt.plot(data[:,0]/3.16,data[:,1]) In the case of the monolayer (figure 2.a) I obtain the same results. thanks