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)
    print(dim)
    zero_array = np.zeros((dim, dim), dtype=float)
    def selfenergy_func(energy, args=()):
        return zero_array

    vlead = kwant.builder.SelfEnergyLead(selfenergy_func, vlead_interface,())
    sys.leads.append(vlead)


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],\
                           check_hermiticity=False,params=par).data
 G21=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[3],out_leads=[2],\
                           check_hermiticity=False,params=par).data


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?
Thanks,
Sayandip