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

Reply via email to