Dear Sayandip,

The Green's function you obtain is a submatrix of the total retarded Green's function. I mean by this that you obtain the Green's function matrix between the two sites that are at L and L-1 as you wrote your code.
I suppose here that the fictitious lead you added on sites L and L-1 are ranked 2 and 3

I hope this helps,

On Tue, Nov 9, 2021 at 7:21 PM Sayandip Dhara <> wrote:
Hi Abbout,
Thanks for pointing out the mistakes that I was making. Now, I want to implement a formula like the one mentioned in the thread before to calculate the equilibrium dc current for a Josephson junction:
I = 2*(KbT/hbar)/sum_{w_p} [trace(H21*G^R_{12}) - trace(H12*G^R_{21})
where w_p are the Matsubara frequencies.
I wrote the following to get the green's functions:
def mount_vlead(sys, vlead_interface, norb):
    dim = norb*len(vlead_interface)
    zero_array = np.zeros((dim, dim), dtype=float)
    def selfenergy_func(energy, args=()):
        return zero_array

    vlead = kwant.builder.SelfEnergyLead(selfenergy_func, vlead_interface,())

lead2 = mount_vlead(syst,[lat(L-1)], 2)
lead3 = mount_vlead(syst,[lat(L)], 2)

syst =syst.finalized()

 G12=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[2],out_leads=[3],\
 G21=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[3],out_leads=[2],\

H12=syst.hamiltonian_submatrix(to_sites=[L-1], from_sites=[L],params=par)
H21=syst.hamiltonian_submatrix(to_sites=[L], from_sites=[L-1],params=par)
So my question is now, that if the Green's functions that I am using above are truly the retarded Green's function of the system?

Abbout Adel