Dear Wilson, Check the code below. (this is not the efficient way to calculate the transport coefficients) I hope this helps, Adel ########################################## import kwant lat = kwant.lattice.square() def make_system(r=5, 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 -4 < pos[0] < 4 lead = kwant.Builder( kwant.TranslationalSymmetry((0, 1))) lead[lat.shape(lead_shape, (0, 0))] = 0 lead[lat.neighbors()] = -t syst.attach_lead(lead) syst.attach_lead(lead.reversed()) return syst syst = make_system() fsyst = syst.finalized() kwant.plot(fsyst), #Get H: H = fsyst.hamiltonian_submatrix(sparse=False) #get the green's function matrix at energy E=-1.2 energy=-1.2 greens_function=kwant.greens_function(fsyst,energy=energy).data print('The shape of the greens function matrix is:' ,greens_function.shape) # index0 = fsyst.lead_interfaces[0] # the index of sites to wich the lead number 0 is connected. index1 = fsyst.lead_interfaces[1] # the index of sites to wich the lead number 1 is connected. print('the total number of sites connected to a lead is: N=', len(index0)+len(index1)) print('system size:', H.shape) print('the greens function is NxN, which is diffferent from the size of the matrix H') def plot_color(site_index): if site_index in index0: return 'g' if site_index in index1: return 'y' else: return 'b' kwant.plot(fsyst,site_color=plot_color), print('If you want the total greens functions that includes the sites that are not connected to a lead') print('you need to construct the matrix manually by adding zero elements at the right place') print('Or you can do it by adding ficticious leads on the other sites. (check the kwant forum)') import numpy as np Sigma=np.zeros(H.shape,complex) # create an initial zero matrix # Fill the element of the matrix with the non-zero part of the self energies (of the leads.) Sigma0 =fsyst.leads[0].selfenergy(energy) Sigma1 =fsyst.leads[1].selfenergy(energy) # we use the fancy indexing in python to call/fill a submatrix Sigma[index0,:][:,index0]=Sigma0 Sigma[index1,:][:,index1]=Sigma1 # the total greens function is therefore # keep in mind that this method is not efficient to calculate the transport coefficients. G=np.linalg.inv(energy -H-Sigma) On Wed, May 29, 2024 at 6:47 AM <wilson2048@outlook.com> wrote:

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.

-- Abbout Adel

Dear Wilson, Sorry, there was a small typo in the previous code: I used (E-H-Sigma) instead of the correct form (E*np.eye(N)-H-Sigma). The code below is now working fine. I hope this helps, Adel ########################################## import kwant import numpy as np lat = kwant.lattice.square() def make_system(r=5, 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 -4 < pos[0] < 4 lead = kwant.Builder( kwant.TranslationalSymmetry((0, 1))) lead[lat.shape(lead_shape, (0, 0))] = 0 lead[lat.neighbors()] = -t syst.attach_lead(lead) syst.attach_lead(lead.reversed()) return syst syst = make_system() fsyst = syst.finalized() kwant.plot(fsyst), #Get H: H = fsyst.hamiltonian_submatrix(sparse=False) #get the green's function matrix at energy E=-1.2 energy=-1.2 greens_function=kwant.greens_function(fsyst,energy=energy).data print('The shape of the greens function matrix is:' ,greens_function.shape) # index0 = fsyst.lead_interfaces[0] # the index of sites to wich the lead number 0 is connected. index1 = fsyst.lead_interfaces[1] # the index of sites to wich the lead number 1 is connected. print('the total number of sites connected to a lead is: N=', len(index0)+len(index1)) print('system size:', H.shape) print('the greens function is NxN, which is diffferent from the size of the matrix H') def plot_color(site_index): if site_index in index0: return 'g' if site_index in index1: return 'y' else: return 'b' kwant.plot(fsyst,site_color=plot_color), print('If you want the total greens functions that includes the sites that are not connected to a lead') print('you need to construct the matrix manually by adding zero elements at the right place') print('Or you can do it by adding ficticious leads on the other sites. (check the kwant forum)') import numpy as np Sigma=np.zeros(H.shape,complex) # create an initial zero matrix # Fill the element of the matrix with the non-zero part of the self energies (of the leads.) Sigma0 =fsyst.leads[0].selfenergy(energy) Sigma1 =fsyst.leads[1].selfenergy(energy) # we use the fancy indexing in python to call/fill a submatrix Sigma[np.ix_(index0, index0)]=Sigma0 Sigma[np.ix_(index1, index1)]=Sigma1 # the total greens function is therefore # keep in mind that this method is not efficient to calculate the transport coefficients. N,_=H.shape G=np.linalg.inv(energy*np.eye(N)-H-Sigma) index=list(index0)+list(index1) print (np.allclose(greens_function,G[np.ix_(index, index)] )) On Thu, May 30, 2024 at 9:02 AM <wilson2048@outlook.com> wrote:

Dear Adel:

Thanks for your reply and the code. I think your code is basically the same as mine, except that I used only one lead and energy=0 to illustrate the problem. My issue is that the corresponding Green's function calculated by inverse is not equal to the one given by kwant. I indeed compared the selected sites, in your language, the surface Green's function.

I hope you can point out why the code of my last line give no zero matrix. Thank you very much.

Bests, Wilson

-- Abbout Adel