Hello, 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 #hv.extension('bokeh') 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. #### syst.attach_lead(lead0) syst.attach_lead(lead1) return syst,lat ---------------------------------------------------------------------------------------------------------------------------- mu=1 t=1.0 delta=0.5 V_N=1 phi=0 L=10 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,()) sys.leads.append(vlead) 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],\ check_hermiticity=False,params=par).data G21=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[L],out_leads=[L-1],\ check_hermiticity=False,params=par).data ---------------------------------------------------------------------------------------------------------------------------------- The error message I am getting is " ValueError: Lead 2 is attached to a site that does not belong to the scattering region: 9" It would be great if you can suggest to me where I am going wrong. Thanks, Sayandip