Dear all,
thanks in advance for the support.
As I already had occasion to explain to Xavier, my problem concerns a
ring geometry where changing Hamiltonian values on non-existing sites
inside the ring (r<r1) actually changes the conductance behavior, while
it should not.
More precisely, I am building a ring, assigning to the ring sites a
certain value given by the/onsite/ function, which however in its
definition contains a conditional /if-/instruction involving sites
inside the ring.
I was expecting that the /if-/condition is then never satisfied and thus
no change is made, but the conductance curve changes if I change the
values of, say, the chemical potential (on-site energy) inside the ring,
where no sites exist.
However, If I do the /same/ changes /outside/ the ring (r>r2), no change
occurs in the conductance, as it should.
Of course one hypothesis is that I am in fact changing also Hamiltonian
elements of the ring, but I don't see why.
As suggested by Xavier, I will try to check that by using plotter.
I attach the piece of code.
I put a comment on the line where I write the conditional instruction
for mu.
If you try to run the code for different values of mu_inside, and/or
different physical sub-regions inside the ring, you should see different
conductance behavior.
I am using python 2.7.5 on a Mac OS X 10.6.8 and the same holds true for
a Mac OS X 10.8.5.
Thanks a lot in advance for any hint on how to understand this misbehavior.
Diego
==============================================================================
import kwant
import scipy, scipy.io
from numpy import *
from pylab import *
import sys
from numpy import mod
from math import pi,atan2
import matplotlib as mpl
from matplotlib import pyplot
import random; import pickle
params = {'backend': 'pdf',
'ytick.direction': 'out',
'text.usetex': True,
'font.size':16}
rcParams.update(params)
matplotlib.rcParams.update({'font.size': 20})
im=0+1j
import tinyarray
import scipy.sparse.linalg
class SimpleNamespace(object):
"""A simple container for parameters."""
def __init__(self, **kwargs):
self.__dict__.update(kwargs)
s_0 = identity(2)
s_z = array([[1, 0], [0, -1]])
s_x = array([[0, 1], [1, 0]])
s_y = array([[0, -1j], [1j, 0]])
#=================================================================================
def make_system(a=1, W=10, r1=10, r2=20):
lat = kwant.lattice.square(a)
sys = kwant.Builder()
def ring(pos):
(x, y) = pos
rsq = x ** 2 + y ** 2
return (r1 ** 2 < rsq < r2 ** 2)
def mu(site, p):
(x, y) = site.pos
if abs(x) < r1 and abs(y) < r1:
return p.mu_inside # this is the critical line: changing
mu values inside the ring should not change the ring conductance
else:
return p.mu_ring
def onsite(site, p):
return ( 4.0 * p.t - mu(site, p)) * s_0
def hopping(site0, site1, p):
return -p.t * s_0
#
===============================================================================
sys[lat.shape(ring, (0, r1 + 1))] = onsite
sys[lat.neighbors()] = hopping
#
===================================================================================
sym_lead = kwant.TranslationalSymmetry((-a, 0))
lead = kwant.Builder(sym_lead)
def lead_shape(pos):
(x, y) = pos
return (-W / 2 < y < W / 2)
lead[lat.shape(lead_shape, (0, 0))] = onsite
lead[lat.neighbors()] = hopping
sys.attach_lead(lead)
sys.attach_lead(lead.reversed())
return sys
def plot_conductance(sys, energies, p):
data = []
for energy in energies:
smatrix = kwant.smatrix(sys, energy, args=[p])
data.append(smatrix.transmission(1, 0))
axes([0.12, 0.12, 0.8, 0.8])
hold(True)
pyplot.plot(energies, data, 'bo-', ms=2.8,
linewidth=0.5,markeredgewidth=0.03)
pyplot.axis([0, 3.0, 0, 10.1])
pyplot.xlabel(r'$eV[t]$')
pyplot.ylabel(r'$G[e^2/h]$')
savefig('Ring_G(e).pdf')
def main():
sys = make_system()
kwant.plot(sys)
sys = sys.finalized()
params = SimpleNamespace(t=1.0, mu_ring=0.0, mu_inside=-1.62)
pyplot.figure()
plot_conductance(sys, energies=[0.01 * i for i in xrange(300)],
p=params)
if __name__ == '__main__':
main()
==================================================================================

Dear all,
I have two simple questions concerning electronic transport with spins with
no interaction:
If I take the program spin_orbit.py from the tutorial and call the system
sys=make_system ( alpha=0, ez=0, W=2, L=2)
# no spin orbit interaction
print matrix.round(kwant.smatrix(sys,energy=1.7).data , 2 )
I choose the energy E=1.7 because one mode is closed.
the result is:
[[-0.00-0.j 0.00-0.j 0.41+0.43j -0.74-0.33j]
[ 0.00-0.j -0.00-0.j -0.80+0.08j -0.52+0.28j]
[ 0.41+0.43j -0.74-0.33j -0.00-0.j 0.00-0.j ]
[-0.80+0.08j -0.52+0.28j 0.00-0.j -0.00-0.j ]]
What surprises me is that there is scattering between spins knowing that
the Hamiltonians for the two spins are separated without interaction !
This seems to appear only when one mode is closed.
To understand better and check this result, I rewrote the program using
two lattices (as for electrons and holes). the program is attached with
the Email.
The result for the scattering matrix is
[[ 0.00-0.j -0.16+0.99j 0.00+0.j 0.00+0.j ]
[-0.16+0.99j 0.00+0.j 0.00+0.j 0.00+0.j ]
[ 0.00+0.j 0.00+0.j 0.00+0.j -0.16+0.99j]
[ 0.00+0.j 0.00+0.j -0.16+0.99j 0.00+0.j ]]
The result is now better because there is no scattering between spins.
My Question is :what am I missing?
The second question concerns the position of the coefficients in the
scattering matrix.
We have the Hamiltonian [image: H=h \otimes 1_{2\times2}] where
[image:
\otimes] is the Kronecker product and [image: h] the spinlessHamiltonian.
The scattering matrix is given by :
[image: S=-1+2\pi i W^\dagger \frac{1}{E-H-\Sigma} W]
[image: W] is the coupling matrix and [image: \Sigma] the self energy
If we look at the simple system [image: ( L=2,W=2)] all the matrices are
square and we can easily proove, just by the properties of the Kronecker
poduct
that the scattering matrix can be written as:
[image: S=s\otimes 1_{2\times2}]
Where [image: s] is the scattering matrix without spin.
[image: s= \left(\begin{array} {c c} r & t \\ t^\prime & r^\prime
\end{array}\right)]
[image: t] is the transmission matrix.
So I am expecting a scattering matrix with the followin form:
[image: S= \left(\begin{array} {c c c c} r &0 & t & 0 \\ 0 &r & 0 & t \\
t^\prime &0 & r^\prime & 0 \\ 0& t^\prime & 0 & r^\prime
\end{array}\right)]
This is consistent with the result with two lattices but still the position
of the coefficients are different from the one obtained
from the Hamiltonian base. So I suggest that you are using different base.
The question I ask is: in order to know the position of coefficients
corresponding to different modes or different spins, which formula
do you use to obtain the scattering matrix. is it a matrix inversion or
some recursive method ?
Thank you in advance.
Regards,
A. Abbout

Hi again,
I just realized that my example was kind of broken as I defined a
a variable `mode` and then promptly overwrote this with something
else. The examples should have read as follows::
lead = 1 # whatever
mode = 0 # whatever
smatrix = kwant.smatrix(sys, energy)
modes = smatrix.lead_info
sites = fsys.leads[lead].sites
psi = modes[lead].wave_functions[:, mode]
psi = dict(zip(sites, psi))
and::
sites = fsys.leads[lead].sites
psi = modes[lead].wave_functions[:, mode]
psi_up = dict(zip(sites, psi[::2])) # even elements of psi (up)
psi_down = dict(zip(sites, psi[1::2])) # odd elements of psi (down)
Regards,
Joe

Hi,
An addition to Joseph’s reply:
If all you want to do is to save the wave function, you will probably
not want to build a dictionary, but simply loop over the wave function
array and the list of sites while writing both to a file. For example
like this (with spin):
sites = fsys.leads[lead].sites
psi = modes[lead].wave_functions[:, mode_nr]
for site, value in zip(sites, psi.reshape(-1, 2)):
tag = site.tag
print tag[0], tag[1], value[0], value[1]
Note how I used NumPy’s “reshape” to be able to iterate over the spins.
Also, bear in mind that the coordinates of the sites that make up a lead
are most likely not the ones that you have set, but equivalent ones
(under that leads’ symmetry).
Best,
Christoph

Dear all,
I would like to ask about the wave function of the propagating modes because I don't know how should i write it into the file. I do that using the code
smatrix = kwant.smatrix(sys, energy)
mode= smatrix.lead_info
and now mode[1].wave_functions[i][0] gives as the wave function of the propagating mode corresponding to the 0 orbital in lead 1. I would like to write this wave function into the file but a don't know how is the correspondence between position "i" and the real position in the lead. It is important for me because I need to know how is the average spin of each propagating mode.
Thank you in advance
Best Regards
Pawel

Dear all,
I am trying square lattice with spin-orbit coupling and Zeeman splititng
just like the example of Tutorial 2.3.1. Now I want to get the
spin-dependent conductance：Guu,Gud,Gdu,and Gdd, where "u" indicates up
spin, d indicates down spin. We can have the scattering matrix through
Smatrix=kwant.smatrix(sys, energy), but how the spin is positioned in the
scattering matrix? I have checked some simple cases, it seems the Smatrix
is organized like this:
| r t |
S= | t' r' |
but how the spin is included?
Thanks in advance.
Qingtiian Zhang
The code is:
import kwant
from matplotlib import pyplot
from numpy import matrix
import tinyarray
sigma_0 = tinyarray.array([[1, 0], [0, 1]])
sigma_x = tinyarray.array([[0, 1], [1, 0]])
sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
sigma_z = tinyarray.array([[1, 0], [0, -1]])
sys=kwant.Builder()
lat=kwant.lattice.square()
a=1
t=1.0
alpha=0.5
e_z=0.08
W=2
L=2
sys[(lat(x, y) for x in range(L) for y in range(W))] = \
4 * t * sigma_0 + e_z * sigma_z
# hoppings in x-direction
sys[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
-t * sigma_0 - 1j * alpha * sigma_y
# hoppings in y-directions
sys[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
-t * sigma_0 + 1j * alpha * sigma_x
#### Define the left lead. ####
lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
lead[(lat(0, j) for j in xrange(W))] = 4 * t * sigma_0 + e_z * sigma_z
# hoppings in x-direction
lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
-t * sigma_0 - 1j * alpha * sigma_y
# hoppings in y-directions
lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
-t * sigma_0 + 1j * alpha * sigma_x
#### Attach the leads and return the finalized system. ####
sys.attach_lead(lead)
sys.attach_lead(lead.reversed())
kwant.plot(sys)
sys=sys.finalized()
print "Hamiltonian="
print sys.hamiltonian_submatrix()
Smatrix=kwant.smatrix(sys, energy=3.)
print "The scattering matrix="
print matrix.round((Smatrix.data),2)
print "T=",(Smatrix.transmission(1, 0))

Hi there,
I'm trying to simulate a scattering region confined between two
reservoirs(leads) for that purpose I defines a wire with linearly increasing
onsite potentials. Now I want to find the local current density at each site
on the wire but I don't know how to do it. Can anyone guide me through this?
:)
def make_system(a=1, t=1.0, W = 60, L = 100, pot=-0.2):
lat = kwant.lattice.square(a)
sys = kwant.Builder()
def hop(site1, site2, B=0):
# The magnetic field is controlled by the parameter B
x1, y1 = site1.pos
x2, y2 = site2.pos
return -t * exp(0.5j * B * (y1+y2)* (x1-x2))#exp(0.5j * B * (y1-x1))
def potential_prof(i,L,pot):
#potential profile linearly determined here
return (L-i)*pot/L
# Define the scattering region
for i in xrange(L):
for j in xrange(W):
# On-site Hamiltonian
sys[lat(i, j)] = 4 * t+potential_prof(i,L,pot)
sys[kwant.builder.HoppingKind((1, 0), lat, lat)] = hop
sys[kwant.builder.HoppingKind((0, 1), lat, lat)] = hop
sym_left_lead = kwant.TranslationalSymmetry((-a, 0))
left_lead = kwant.Builder(sym_left_lead)
for j in xrange(W):
left_lead[lat(0, j)] = 4 * t+pot
if j > 0:
left_lead[lat(0, j), lat(0, j - 1)] = -t
left_lead[lat(1, j), lat(0, j)] = -t
sys.attach_lead(left_lead)
sym_right_lead = kwant.TranslationalSymmetry((a, 0))
right_lead = kwant.Builder(sym_right_lead)
for j in xrange(W):
right_lead[lat(0, j)] = 4 * t
if j > 0:
right_lead[lat(0, j), lat(0, j - 1)] = -t
right_lead[lat(1, j), lat(0, j)] = -t
sys.attach_lead(right_lead)
# Modify the scattering region
del sys[lat(85, 2)]
del sys[lat(87, 2)]
del sys[lat(86, 3)]
del sys[lat(25, 37)]
del sys[lat(27, 37)]
del sys[lat(26, 38)]
sys = sys.finalized()
return sys