Dear Adel: Thanks for your reply and the code. I have read it through (also the linked thread). I am aware that the green's function is not the total one.
What I am doing is obtain the total one by inverse the whole Hamiltonian matrix plus the self-energy. Then select the part that connected to the leads, and compare it with the one kwant calculated. I found they don't match. I have added the some comments on my code to clarify my question ``` lat = kwant.lattice.square() def make_system(r=20, t=-1): def circle(pos): x, y = pos return x**2 + y**2 < r**2 syst = kwant.Builder() syst[lat.shape(circle, (0, 0))] = 0 syst[lat.neighbors()] = t syst.eradicate_dangling() def lead_shape(pos): return -15 < pos[0] < 15 lead = kwant.Builder( kwant.TranslationalSymmetry((0, 1))) lead[lat.shape(lead_shape, (0, 0))] = 0 lead[lat.neighbors()] = -t syst.attach_lead(lead) return syst syst = make_system() fsyst = syst.finalized() kwant.plot(fsyst) greens_function=kwant.greens_function(fsyst) Hc = fsyst.hamiltonian_submatrix(sparse=False) # the index of the sites IN SYSTEM which connected to lead index = fsyst.lead_interfaces[0] # add the self energy to the original matrix Hc[np.ix_(index, index)] += greens_function.lead_info[0] # inverse it, to get the total Green's function of the whole system Gr = np.linalg.inv(-Hc) # Should be a zero matrix! The first term is the green's function of the selected sites. Gr[np.ix_(index, index)]-greens_function.submatrix(0, 0) ``` I am also very interested in your last comment, about how to get the whole green's function of the system (include the bulk sites). I hope you can provide me some links, thank you very much.