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
