Dear kwant users I'm trying to implement the interband-kubo formula for dc-conductivity.
The underlying system is a square lattice with attached leads in z-direction and perpendicular magnetic field. My problem is that I couldn't manage to calculate the expectation value for the current operator for different bands. Could anyone give me a hint what is the right way of doing it? Best regards Simon def fermi(energy, Efermi): f = np.real(1 / (exp((energy - Efermi) / temp) + 1)) return f def make_system(a=1, t=1.0, w=15, d=60): lat = kwant.lattice.square(a, norbs=1) syst = kwant.Builder() def hopping(site_i, site_j, phi): xi, yi = site_i.pos xj, yj = site_j.pos if yi == yj: return -t if xi == xj and yi != yj: return (-0.01 * t) * exp(-0.5j * phi * (xi + xj)) syst[(lat(x, y) for x in range(w) for y in range(d))] = 4 * t syst[lat.neighbors()] = hopping lead = kwant.Builder(kwant.TranslationalSymmetry((0, -a))) lead[(lat(j, 0) for j in range(w))] = 4 * t lead[lat.neighbors()] = hopping syst.attach_lead(lead) syst.attach_lead(lead.reversed()) return syst syst = make_system() syst = syst.finalized() def Interband_Kubo(syst): w=15 d=60 Gs_inter = [] count = 0 phis = np.linspace(0, (2 * pi / 15) * 5, 100) for phi in phis: count += 1 bands = kwant.physics.Bands(syst.leads[1], params=dict(phi=phi)) momenta = np.linspace(-pi, pi, 300) energydisp = [bands(k) for k in momenta] elist = [] for m in range(0, len(energydisp), 1): for n in range(0, len(energydisp[m]), 1): energy = np.real(energydisp[m][n]) elist.append(energy) selist = np.sort(elist) EF = selist[int(len(elist) / 2)] condsum = [] for l in range(0, len(energydisp), 1): for s in range(0, len(energydisp[l]), 1): count = 0 for q in range(0, len(energydisp[l]), 1): count += 1 if s != q: #print(count) energy_s = np.real(energydisp[l][s]) energy_q = np.real(energydisp[l][q]) if EF -40*temp <= energy_s <= EF + 40*temp: if EF - 40 * temp <= energy_q <= EF + 40* temp: J_z = kwant.operator.Current(syst) wf_s = kwant.solvers.default.wave_function(syst, energydisp[l][s], params=dict(phi=phi)) wf_q = kwant.solvers.default.wave_function(syst, energydisp[l][q], params=dict(phi=phi)) if len(wf_s(1)) == 1 and len(wf_q(1)) ==1: psi_s = wf_s(1)[0] psi_q = wf_q(1)[0] print(psi_s) J_qs = np.sum(J_z(psi_q, psi_s, params=dict(phi=phi))) J_sq = np.sum(J_z(psi_s, psi_q, params=dict(phi=phi))) else: J_qs = 0 J_sq = 0 cond = 1j * ((fermi(energy_s, EF) - fermi(energy_q, EF))) * ( (J_sq * J_qs) / (energy_s - energy_q + (1 / 15 ** 2) * 1j)) condsum.append(cond) conductance = (np.sum(condsum) / (15*60)) print('conductance', conductance) Gs_inter.append(conductance) return Gs_inter fluxes = np.linspace(0, 5, 100) plt.plot(fluxes, Interband_Kubo(syst), label = "inter") plt.grid(True) plt.title(" Conductance") plt.legend() plt.ylabel("conductance") plt.xlabel(r'$\phi [2\pi / L]$') plt.show()