I am trying to calculate the Green's function for each layer/slice of a ribbon/wire. I can find the self energy and surface Green's function from which (in principle, I have not tried yet) it is possible find the Green's function for successive layer using RGF algorithm as discussed in this thread (http://thread.gmane.org/gmane.comp.science.kwant.user/647).
From the GreensFunction documentation, it looks like there is a way to specify the slice/block of the system and directly find the Green's function, but I don't find any example of that. Is there any way to do that ? ( http://kwantproject.org/doc/1/reference/generated/kwant.solvers.common.Gree... )
For example consider this graphene ribbon. I added a random impurity to break the translational symmetry so that I can distinguish each slice. I can find the Green's function for the sites which are connected to leads only. How do I calculate the Green's function for each slice? Best, Sumit  from __future__ import division from math import sqrt from matplotlib import pyplot import kwant import numpy as np sin_30, cos_30 = (1 / 2, sqrt(3) / 2) graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)], [(0, 0), (0, 1 / sqrt(3))]) a, b = graphene.sublattices L=3;W=3; def box(pos): #scattering region x, y = pos return L<x<L and W<y<W def lshape(pos): #lead x,y = pos return W<y<W def onsite(site): return kwant.digest.uniform(repr(site)) hoppings = (((0, 0), a, b), ((0, 1), a, b), ((1, 1), a, b)) def scater(): sys = kwant.Builder() sys[graphene.shape(box, (0, 0))] = onsite sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = 1 lead = kwant.Builder(kwant.TranslationalSymmetry(graphene.vec((1, 0)))) lead[graphene.shape(lshape, (0, 0))] = 0 lead[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = 1 sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys sys=scater() #def family_colors(site): # return 0 if site.family == a else 1 #kwant.plot(sys, site_color=family_colors, site_lw=0.1, colorbar=False) sys = sys.finalized() # Self energy flead0 = sys.leads[0] flead1 = sys.leads[1] s_en = flead0.selfenergy(0.5) #Green's function g = kwant.greens_function(sys,0.5) gg = g.submatrix(1,0)  Sumit Ghosh Spintronics Theory Group PSE Division KAUST   This message and its contents, including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.
Dear Sumit, The easiest way to calculate Green’s function on arbitrary set of sites in your system is to specify additional virtual lead with zero self energy with following builder class: http://kwantproject.org/doc/1/reference/generated/kwant.builder.SelfEnergyL... Here is an example function, that can do it: def mount_vlead(sys, vlead_interface, norb): """Mounts virtual lead to interfaces provided. :sys: kwant.builder.Builder An unfinalized system to mount leads :vlead_interface: sequence of kwant.builder.Site Interface of lead :norb: integer Number of orbitals in system hamiltonian. """ dim = len(vlead_interface)*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) Then you can calculate Green’s gunction on interface of this lead. Basis in that case is specified on order of sites in vlead_interface. If I remember correctly, calculation time in that case scales as n*N*log(N), where N is number of sites in system, n is number of sites to calculate Green’s function in. With best regards, Viacheslav Ostroukh InstituutLorentz – Niels Bohrweg 2 – Room 259 – 2333 CA Leiden ostroukh@ilorentz.org slava@ostroukh.me +31 6 444 968 12 +38 099 721 76 06 From: Sumit Ghosh Sent: Tuesday, January 12, 2016 11:44 To: kwantdiscuss@kwantproject.org Subject: [Kwant] Green's function for each block I am trying to calculate the Green's function for each layer/slice of a ribbon/wire. I can find the self energy and surface Green's function from which (in principle, I have not tried yet) it is possible find the Green's function for successive layer using RGF algorithm as discussed in this thread (http://thread.gmane.org/gmane.comp.science.kwant.user/647).
From the GreensFunction documentation, it looks like there is a way to specify the slice/block of the system and directly find the Green's function, but I don't find any example of that. Is there any way to do that ? (http://kwantproject.org/doc/1/reference/generated/kwant.solvers.common.Gree...)
For example consider this graphene ribbon. I added a random impurity to break the translational symmetry so that I can distinguish each slice. I can find the Green's function for the sites which are connected to leads only. How do I calculate the Green's function for each slice? Best, Sumit  from __future__ import division from math import sqrt from matplotlib import pyplot import kwant import numpy as np sin_30, cos_30 = (1 / 2, sqrt(3) / 2) graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)], [(0, 0), (0, 1 / sqrt(3))]) a, b = graphene.sublattices L=3;W=3; def box(pos): #scattering region x, y = pos return L<x<L and W<y<W def lshape(pos): #lead x,y = pos return W<y<W def onsite(site): return kwant.digest.uniform(repr(site)) hoppings = (((0, 0), a, b), ((0, 1), a, b), ((1, 1), a, b)) def scater(): sys = kwant.Builder() sys[graphene.shape(box, (0, 0))] = onsite sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = 1 lead = kwant.Builder(kwant.TranslationalSymmetry(graphene.vec((1, 0)))) lead[graphene.shape(lshape, (0, 0))] = 0 lead[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = 1 sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys sys=scater() #def family_colors(site): # return 0 if site.family == a else 1 #kwant.plot(sys, site_color=family_colors, site_lw=0.1, colorbar=False) sys = sys.finalized() # Self energy flead0 = sys.leads[0] flead1 = sys.leads[1] s_en = flead0.selfenergy(0.5) #Green's function g = kwant.greens_function(sys,0.5) gg = g.submatrix(1,0)  Sumit Ghosh Spintronics Theory Group PSE Division KAUST This message and its contents, including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.
Thanks Viacheslav. Now I am stuck with the vlead_interface part. According to the documentation it is a sequence of Site instances, which as much I understand would be a list of the position and SiteFamily, but I don't know how to define them. For example, in the earlier graphene nanoribbon case, let say I want to choose this slice x==1, y in range(W,W) or I can define a region like def vshape(pos): #vlead x,y = pos return x==1 and W<y<W How do I define vlead_interface here? Best regards, Sumit On 12 January 2016 at 13:37, Viacheslav Ostroukh < ostroukh@lorentz.leidenuniv.nl> wrote:
Dear Sumit,
The easiest way to calculate Green’s function on arbitrary set of sites in your system is to specify additional virtual lead with zero self energy with following builder class:
http://kwantproject.org/doc/1/reference/generated/kwant.builder.SelfEnergyL...
Here is an example function, that can do it:
def mount_vlead(sys, vlead_interface, norb):
"""Mounts virtual lead to interfaces provided.
:sys: kwant.builder.Builder
An unfinalized system to mount leads
:vlead_interface: sequence of kwant.builder.Site
Interface of lead
:norb: integer
Number of orbitals in system hamiltonian.
"""
dim = len(vlead_interface)*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)
Then you can calculate Green’s gunction on interface of this lead. Basis in that case is specified on order of sites in vlead_interface.
If I remember correctly, calculation time in that case scales as n*N*log(N), where N is number of sites in system, n is number of sites to calculate Green’s function in.
With best regards, Viacheslav Ostroukh
InstituutLorentz – Niels Bohrweg 2 – Room 259 – 2333 CA Leiden
ostroukh@ilorentz.org slava@ostroukh.me
+31 6 444 968 12 +38 099 721 76 06
*From: *Sumit Ghosh <sumit.ghosh@kaust.edu.sa> *Sent: *Tuesday, January 12, 2016 11:44 *To: *kwantdiscuss@kwantproject.org *Subject: *[Kwant] Green's function for each block
I am trying to calculate the Green's function for each layer/slice of a ribbon/wire. I can find the self energy and surface Green's function from which (in principle, I have not tried yet) it is possible find the Green's function for successive layer using RGF algorithm as discussed in this thread (http://thread.gmane.org/gmane.comp.science.kwant.user/647).
From the GreensFunction documentation, it looks like there is a way to specify the slice/block of the system and directly find the Green's function, but I don't find any example of that. Is there any way to do that ?
( http://kwantproject.org/doc/1/reference/generated/kwant.solvers.common.Gree... )
For example consider this graphene ribbon. I added a random impurity to break the translational symmetry so that I can distinguish each slice. I can find the Green's function for the sites which are connected to leads only. How do I calculate the Green's function for each slice?
Best,
Sumit

from __future__ import division
from math import sqrt
from matplotlib import pyplot
import kwant
import numpy as np
sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
[(0, 0), (0, 1 / sqrt(3))])
a, b = graphene.sublattices
L=3;W=3;
def box(pos): #scattering region
x, y = pos
return L<x<L and W<y<W
def lshape(pos): #lead
x,y = pos
return W<y<W
def onsite(site):
return kwant.digest.uniform(repr(site))
hoppings = (((0, 0), a, b), ((0, 1), a, b), ((1, 1), a, b))
def scater():
sys = kwant.Builder()
sys[graphene.shape(box, (0, 0))] = onsite
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = 1
lead = kwant.Builder(kwant.TranslationalSymmetry(graphene.vec((1, 0))))
lead[graphene.shape(lshape, (0, 0))] = 0
lead[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = 1
sys.attach_lead(lead)
sys.attach_lead(lead.reversed())
return sys
sys=scater()
#def family_colors(site):
# return 0 if site.family == a else 1
#kwant.plot(sys, site_color=family_colors, site_lw=0.1, colorbar=False)
sys = sys.finalized()
# Self energy
flead0 = sys.leads[0]
flead1 = sys.leads[1]
s_en = flead0.selfenergy(0.5)
#Green's function
g = kwant.greens_function(sys,0.5)
gg = g.submatrix(1,0)

Sumit Ghosh
Spintronics Theory Group
PSE Division
KAUST
This message and its contents, including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.
 Sumit Ghosh Spintronics Theory Group PSE Division KAUST   This message and its contents, including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.
Hi,
Now I am stuck with the vlead_interface part. According to the documentation it is a sequence of Site instances, which as much I understand would be a list of the position and SiteFamily, but I don't know how to define them.
`Site` objects are returned by lattice objects in Kwant. whenever you do: lat = kwant.lattice.square() site = lat(0, 1) then `site` is an instance of a `Site` with `lat` as its family and a tag of `(0, 1)`. You need to provide the interface sites so that Kwant knows how you have ordered the degrees of freedom in your self energy matrix. Given that you are creating a virtual lead and want to get the Green's function everywhere in the system, you could do something like: syst = kwant.Builder() ... # build the system greens_function_sites = syst.sites() mount_vlead(syst, greens_function_sites, 1) fsyst = syst.finalized() The only caveat is that when you obtain the Green's function at the interface with the virtual lead (i.e. on all the sites) the degrees of freedom will be ordered like the `greens_function_sites` sequence, and *not* in the same way as `fsyst.sites`. Hope that helps, Joe
Thanks Joe for your suggestion. I am able to calculate the Green's function at any arbitrary position, but I can't figure out the sequence. How can I find the order of site and orbital corresponding to each element of the Green's function matrix? For example,  def scater() make system and attach lead, not finalized # mounting a virtual lead #Viacheslav def mount_vlead(sys, vlead_interface, norb): dim = len(vlead_interface)*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) # arbitrary virtual lead #Joe def vshape(pos): x, y = pos return 0.5<=x<2.5 and 0.5<=y<=W vlead = kwant.Builder(kwant.TranslationalSymmetry(graphene.vec((1, 0)))) vlead[graphene.shape(vshape, (1.5,0.5))] = 0 vlead[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = 1 gf_sites = vlead.sites() syst = scater() mount_vlead(syst, gf_sites, 1) fsyst = syst.finalized()  When I simply print gf_sites it shows all the information as dict_keys(), I can not extract the data from that. I tried to use kwant.plotter.sys_leads_sites(fsyst, num_lead_cells=1) as described in the documentation, but can not make anything there as well. What shall I do? There is another thing I would like to clarify. Since my objective is to find the Green's function at each lattice points, would it be correct to reduce vshape such that it contains only one lattice point? By doing so I may destroy the geometry of the virtual lead, but we are defining the self energy to be a null matrix. So that should not affect the Green's function. Thanks and best regards, Sumit On 12 January 2016 at 19:46, Joseph Weston <joseph.weston08@gmail.com> wrote:
Hi,
Now I am stuck with the vlead_interface part. According to the documentation it is a sequence of Site instances, which as much I understand would be a list of the position and SiteFamily, but I don't know how to define them.
`Site` objects are returned by lattice objects in Kwant. whenever you do:
lat = kwant.lattice.square() site = lat(0, 1)
then `site` is an instance of a `Site` with `lat` as its family and a tag of `(0, 1)`. You need to provide the interface sites so that Kwant knows how you have ordered the degrees of freedom in your self energy matrix. Given that you are creating a virtual lead and want to get the Green's function everywhere in the system, you could do something like:
syst = kwant.Builder() ... # build the system greens_function_sites = syst.sites() mount_vlead(syst, greens_function_sites, 1) fsyst = syst.finalized()
The only caveat is that when you obtain the Green's function at the interface with the virtual lead (i.e. on all the sites) the degrees of freedom will be ordered like the `greens_function_sites` sequence, and *not* in the same way as `fsyst.sites`.
Hope that helps,
Joe
 Sumit Ghosh Spintronics Theory Group PSE Division KAUST   This message and its contents, including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.
participants (3)

Joseph Weston

Sumit Ghosh

Viacheslav Ostroukh