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()
