Creating virtual self-energy leads to extract Green's function between arbitrary sites:
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
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],\ check_hermiticity=False,params=par).data 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],\ check_hermiticity=False,params=par).data You can check the previous post for more details [1] I hope this helps, Adel [1] https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01356.html On Thu, Oct 21, 2021 at 7:20 PM Sayandip Dhara <sayandip@knights.ucf.edu> wrote:
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
-- Abbout Adel
Hi Abbout, Thanks for pointing out the mistakes that I was making. Now, I want to implement a formula like the one mentioned in the thread before to calculate the equilibrium dc current for a Josephson junction: I = 2*(KbT/hbar)/sum_{w_p} [trace(H21*G^R_{12}) - trace(H12*G^R_{21}) where w_p are the Matsubara frequencies. I wrote the following to get the green's functions: _______________________________________________________________ def mount_vlead(sys, vlead_interface, norb): dim = norb*len(vlead_interface) print(dim) 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) lead2 = mount_vlead(syst,[lat(L-1)], 2) lead3 = mount_vlead(syst,[lat(L)], 2) syst =syst.finalized() G12=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[2],out_leads=[3],\ check_hermiticity=False,params=par).data G21=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[3],out_leads=[2],\ check_hermiticity=False,params=par).data H12=syst.hamiltonian_submatrix(to_sites=[L-1], from_sites=[L],params=par) H21=syst.hamiltonian_submatrix(to_sites=[L], from_sites=[L-1],params=par) So my question is now, that if the Green's functions that I am using above are truly the retarded Green's function of the system? Thanks, Sayandip
Dear Sayandip, The Green's function you obtain is a submatrix of the total retarded Green's function. I mean by this that you obtain the Green's function matrix between the two sites that are at L and L-1 as you wrote your code. I suppose here that the fictitious lead you added on sites L and L-1 are ranked 2 and 3 I hope this helps, Adel On Tue, Nov 9, 2021 at 7:21 PM Sayandip Dhara <sayandip@knights.ucf.edu> wrote:
Hi Abbout, Thanks for pointing out the mistakes that I was making. Now, I want to implement a formula like the one mentioned in the thread before to calculate the equilibrium dc current for a Josephson junction: I = 2*(KbT/hbar)/sum_{w_p} [trace(H21*G^R_{12}) - trace(H12*G^R_{21}) where w_p are the Matsubara frequencies. I wrote the following to get the green's functions: _______________________________________________________________ def mount_vlead(sys, vlead_interface, norb): dim = norb*len(vlead_interface) print(dim) 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)
lead2 = mount_vlead(syst,[lat(L-1)], 2) lead3 = mount_vlead(syst,[lat(L)], 2)
syst =syst.finalized()
G12=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[2],out_leads=[3],\ check_hermiticity=False,params=par).data G21=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[3],out_leads=[2],\ check_hermiticity=False,params=par).data
H12=syst.hamiltonian_submatrix(to_sites=[L-1], from_sites=[L],params=par) H21=syst.hamiltonian_submatrix(to_sites=[L], from_sites=[L-1],params=par) So my question is now, that if the Green's functions that I am using above are truly the retarded Green's function of the system? Thanks, Sayandip
-- Abbout Adel
participants (2)
-
Abbout Adel -
Sayandip Dhara