Dear Sayandip,

You have few mistakes in your code:
1) mount_vlead(syst,[L-1], 2). Here you have to provide a list of sites (an interface) on which you want to attach a virtual lead.
   Example: mount_vlead(syst,[lat(1),lat(5),lat(6)], 2)
2) inside the function "mount_vlead" you have to define dim=len(vlead_interface)*norbs.
3) G12=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[L-1],out_leads=[L],\
     Here, you have to provide the rank of the lead following the order you attached them. 0 is the left lead, 1 is the right lead, 2 is the first virtual lead, ....
Example:G12=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[2],out_leads=[3],\

You can check the previous post for more details [1]
I hope this helps,


On Thu, Oct 21, 2021 at 7:20 PM Sayandip Dhara <> wrote:
I am trying to extract Green's function between two arbitrary sites in a system. I have noticed from one of the previous threads that one way to do it is to attach virtual leads with zero self-energy on those specific sites. However, I am having trouble with the suggested function. It is telling me that whichever site I want to attach is not in the scattering region. The following is the code:
import kwant
from types import SimpleNamespace
from copy import copy
import tinyarray
import numpy as np
import csv
import ipywidgets
#import holoviews as hv
from scipy import sparse
from matplotlib import pyplot
from scipy.optimize import minimize
sigma_x = tinyarray.array([[0, 1], [1, 0]])
sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
sigma_z = tinyarray.array([[1, 0], [0, -1]])
sigma_0 = tinyarray.array([[1, 0], [0, 1]])

tau_x = tinyarray.array([[0, 1], [1, 0]])
tau_y = tinyarray.array([[0, -1j], [1j, 0]])
tau_z = tinyarray.array([[1, 0], [0, -1]])
tau_0 = tinyarray.array([[1, 0], [0, 1]])

# mu=0.5,
#                 t=-1.,delta=0.5,V_N=1.25, phi=0, L=3

def make_system_ex3(L=10):
    lat = kwant.lattice.chain(norbs=2)
    syst = kwant.Builder()

    #### Define the scattering region. ####
    def onsite_sc(site,t,mu,delta):
        return (2.*t-mu )*tau_z + delta*tau_x

    def onsite_N(site,t,mu,V_N):
        return (2.*t-mu + V_N)*tau_z

    syst[(lat(x) for x in range(1, L))] = onsite_sc

    #superconducting onsite hamiltonian

    syst[lat(L)] = onsite_N  #normal onsite hamiltonian

    syst[(lat(x) for x in range(L+1,2*L+1))] = onsite_sc #superconducting onsite hamiltonian

    def hop(site1,site2,t,phi):
        return -t*np.matmul(np.exp(1j*phi*tau_z),tau_z)

    syst[(lat(L-1), lat(L))] = hop
    for i in range(1,L-1):
        syst[(lat(i), lat(i+1))] = -t*tau_z
    for i in range(L,2*L):
        syst[(lat(i), lat(i+1))] = -t*tau_z

    #### Define the leads. ####
    sym_left = kwant.TranslationalSymmetry([-1])
    lead0 = kwant.Builder(sym_left)
    lead0[(lat(0))] = (2.*t-mu)*tau_z + delta*tau_x
    lead0[lat.neighbors()] = -t*tau_z

    sym_right = kwant.TranslationalSymmetry([1])
    lead1 = kwant.Builder(sym_right)
    lead1[(lat(2*L+1))] =  (2.*t-mu)*tau_z + delta*tau_x
    lead1[lat.neighbors()] = -t*tau_z

#     #### Attach the leads and return the system. ####

    return syst,lat

par = dict(t=t,mu=mu,delta=delta,phi=phi,V_N=V_N)

syst,lat = make_system_ex3(L=10)

def mount_vlead(sys, vlead_interface, norb):
    dim = norb
    zero_array = np.zeros((dim, dim), dtype=float)
    def selfenergy_func(energy, args=()):
        return zero_array

    vlead = kwant.builder.SelfEnergyLead(selfenergy_func, vlead_interface,())

mount_vlead(syst,[L-1], 2)     # number of orbitals =4. Here we mount lead number 2
mount_vlead(syst,[L], 2)     # number of orbitals =4. Here we mount lead number 3

syst =syst.finalized()

G12=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[L-1],out_leads=[L],\
G21=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[L],out_leads=[L-1],\

The error message I am getting is
" ValueError: Lead 2 is attached to a site that does not belong to the scattering region:
It would be great if you can suggest to me where I am going wrong.

Abbout Adel