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 <wilson2...@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
