Dear All, I am using Kwant to calculate the conductance of a Graphene structure. However, it gives different conductance each time I run the code. I didn't change the code. For the same code, I am getting different results each time. Is there some error in my code or am I using the kwant function incorrectly? Please Help Code ----------- import kwant import matplotlib.pyplot as plt Wnr = 99 lnr = 60 a = 2.46 lat = kwant.lattice.honeycomb(a,norbs=1) a,b = lat.sublattices def make_sys(): def rect(pos): x,y = pos return 0<=x<=lnr and 0<=y<=Wnr def delet(pos): x, y = pos return 0<=x<=35 and 20<=y<=80 model1 = kwant.Builder() model1[lat.shape(rect,(1,1))] = 0 model1[lat.neighbors()] = 2.64 del model1[lat.shape(delet, (11,25))] model1.eradicate_dangling() return model1 def lead_attach(model1): def lead1_shape(pos): x, y = pos return 0<y<=20 def lead2_shape(pos): x, y = pos return 80<=y<=Wnr sym = kwant.TranslationalSymmetry(lat.vec((-1,0))) #from the left sym.add_site_family(lat.sublattices[0], other_vectors=[(-1, 2)]) sym.add_site_family(lat.sublattices[1], other_vectors=[(-1, 2)]) lead = kwant.Builder(sym) lead[lat.shape(lead1_shape, (1,2))]=0 lead[lat.neighbors()] =2.64 model1.attach_lead(lead) ########### lead 2 ################# sym1 = kwant.TranslationalSymmetry(lat.vec((-1,0))) #from the left sym1.add_site_family(lat.sublattices[0], other_vectors=[(-1, 2)]) sym1.add_site_family(lat.sublattices[1], other_vectors=[(-1, 2)]) lead2 = kwant.Builder(sym1) lead2[lat.shape(lead2_shape, (2,82))] = 0 lead2[lat.neighbors()] = 2.64 model1.attach_lead(lead) model1.attach_lead(lead2) kwant.plot(model1) return model1 def dos_eng(model1): rho = kwant.kpm.SpectralDensity(model1) eng1, densities = rho.energies, rho.densities return eng1, densities def calc_conductance(sys, energies): data = [] for energy in energies: smatrix = kwant.smatrix(sys, energy) data.append(smatrix.transmission(1, 0)) plt.figure() plt.plot(energies, data) plt.xlim([-3,3]) plt.ylim([0,4]) plt.xlabel("Energy (eV)",fontweight='bold') plt.ylabel("Conductance $(e^2/h)$",fontweight='bold') return data def main(): sys = make_sys() syst = lead_attach(sys) syst = syst.finalized() eng1, densities = dos_eng(syst) cond = calc_conductance(syst,eng1) if __name__ == '__main__': main()