[Kwant] up/down conductance

[Patrik Arvoy]

Dear kwant users and developers,

I am trying to calculate the up and down conductance in the following
system, separately, to calculate spin polarization.
I followed some of the steps suggested in the mailing list before
(wavefuncation or current operator).
I know for this system I should get the same conductance for up and down
spins.
I first calculate the total conductance
But then
1- The wave function approach does not give the the up/down conductance
similar to total.
2- the current-operator approach gives the error
'kwant.graph.core.EdgeDoesNotExistError'.

I have attached the code in the following. If one runs the code
it plot the system with leads correctly
the first conductance plot is what one expects (if trans:)
then the second plot is incorrect (if strans1:) and then gives the error
(if strans:).

I appreciate any help
Patrik
-----------------------------------------------------------------


from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import *
from matplotlib import rcParams
from numpy import *
from numpy.linalg import *
import pickle
import sys
import os
import string
import heapq
import kwant
import tinyarray
from matplotlib import pyplot

chiral=True
if chiral:
p = pi/5 #phi
t = 0.66 #theta
a = 0.34
x = 1.4
e1 = 0
e2 = 0.3
t2=0.1
t1=-x*t2
t0 = 2
lam=-0.08
t_so1 = 0.01 #spin-orbit coupling param
t_so2 = x*t_so1 #spin-orbit coupling param
tl=tr=0.3
N = 30
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]])
no=2 #number of orbitals
def sigma_v1(ap): # pauli metrix along the vertical axis
value=sigma_z*cos(t)+sin(t)*(sigma_x*sin(ap)-sigma_y*cos(ap))
return value

def sigma_v2(ap): # pauli metrix along the vertical axis
value=sigma_z*cos(t)-sin(t)*(sigma_x*sin(ap)-sigma_y*cos(ap))
return value

def family_color(sites):
return 'black' #if site.family == sites

def hopping_lw(site1, site2):
return 0.08


class Amorphous(kwant.builder.SiteFamily):
def __init__(self, coords):
self.coords = coords
super(Amorphous, self).__init__("amorphous", "",no)

def normalize_tag(self, tag):
try:
tag = int(tag[0])
except:
raise KeyError

if 0 <= tag < len(coords):
return tag
else:
raise KeyError

def pos(self, tag):
return self.coords[tag]

coords=[(0.0000000000, 0.0000000000, 0.0000000000), (-0.1336881039,
0.4114496766, 0.3400000000), (-0.4836881039, 0.6657395614, 0.6800000000),
(-0.9163118961, 0.6657395614, 1.0200000000), (-1.2663118961, 0.4114496766,
1.3600000000), (-1.4000000000, 0.0000000000, 1.7000000000), (-1.2663118961,
-0.4114496766, 2.0400000000), (-0.9163118961, -0.6657395614, 2.3800000000),
(-0.4836881039, -0.6657395614, 2.7200000000), (-0.1336881039,
-0.4114496766, 3.0600000000), (0.0000000000, -0.0000000000, 3.4000000000),
(-0.1336881039, 0.4114496766, 3.7400000000), (-0.4836881039, 0.6657395614,
4.0800000000), (-0.9163118961, 0.6657395614, 4.4200000000), (-1.2663118961,
0.4114496766, 4.7600000000), (-1.4000000000, 0.0000000000, 5.1000000000),
(-1.2663118961, -0.4114496766, 5.4400000000), (-0.9163118961,
-0.6657395614, 5.7800000000), (-0.4836881039, -0.6657395614, 6.1200000000),
(-0.1336881039, -0.4114496766, 6.4600000000), (0.0000000000, -0.0000000000,
6.8000000000), (-0.1336881039, 0.4114496766, 7.1400000000), (-0.4836881039,
0.6657395614, 7.4800000000), (-0.9163118961, 0.6657395614, 7.8200000000),
(-1.2663118961, 0.4114496766, 8.1600000000), (-1.4000000000, 0.0000000000,
8.5000000000), (-1.2663118961, -0.4114496766, 8.8400000000),
(-0.9163118961, -0.6657395614, 9.1800000000), (-0.4836881039,
-0.6657395614, 9.5200000000), (-0.1336881039, -0.4114496766, 9.8600000000),
(-1.4000000000, 0.0000000000, 0.0000000000), (-1.2663118961, -0.4114496766,
0.3400000000), (-0.9163118961, -0.6657395614, 0.6800000000),
(-0.4836881039, -0.6657395614, 1.0200000000), (-0.1336881039,
-0.4114496766, 1.3600000000), (0.0000000000, -0.0000000000, 1.7000000000),
(-0.1336881039, 0.4114496766, 2.0400000000), (-0.4836881039, 0.6657395614,
2.3800000000), (-0.9163118961, 0.6657395614, 2.7200000000), (-1.2663118961,
0.4114496766, 3.0600000000), (-1.4000000000, 0.0000000000, 3.4000000000),
(-1.2663118961, -0.4114496766, 3.7400000000), (-0.9163118961,
-0.6657395614, 4.0800000000), (-0.4836881039, -0.6657395614, 4.4200000000),
(-0.1336881039, -0.4114496766, 4.7600000000), (0.0000000000, -0.0000000000,
5.1000000000), (-0.1336881039, 0.4114496766, 5.4400000000), (-0.4836881039,
0.6657395614, 5.7800000000), (-0.9163118961, 0.6657395614, 6.1200000000),
(-1.2663118961, 0.4114496766, 6.4600000000), (-1.4000000000, 0.0000000000,
6.8000000000), (-1.2663118961, -0.4114496766, 7.1400000000),
(-0.9163118961, -0.6657395614, 7.4800000000), (-0.4836881039,
-0.6657395614, 7.8200000000), (-0.1336881039, -0.4114496766, 8.1600000000),
(0.0000000000, -0.0000000000, 8.5000000000), (-0.1336881039, 0.4114496766,
8.8400000000), (-0.4836881039, 0.6657395614, 9.1800000000), (-0.9163118961,
0.6657395614, 9.5200000000), (-1.2663118961, 0.4114496766, 9.8600000000)]
amorphous_lat = Amorphous(coords)

syst = kwant.Builder()

#adding the onsite and hopping to the system
for i in range(N):
syst[amorphous_lat(i)] = e1*sigma_0
syst[amorphous_lat(N+i)] = e2*sigma_0
syst[amorphous_lat(i), amorphous_lat(N+i)] = lam*sigma_0
if i > 0:
syst[amorphous_lat(i), amorphous_lat(i-1)] = t1*sigma_0 +
1j*t_so1*(sigma_v1(i*p)+sigma_v1((i-1)*p))
syst[amorphous_lat(N+i),amorphous_lat(N+i-1)] = t2*sigma_0 +
1j*t_so2*(sigma_v2(i*p)+sigma_v2((i-1)*p))
# If we want to attach to vertical 1D chains to the system
# we first add a site of the down lead to the scattering region
lat=kwant.lattice.cubic(a, norbs=no)

syst[lat(0, 0, -1)] = e1*sigma_0
syst[amorphous_lat(0), lat(0, 0, -1)] = tl*sigma_0

# We make a regular down lead and attach it to the system
dn_lead = kwant.Builder(kwant.TranslationalSymmetry((0, 0, -a)))
dn_lead[lat(0, 0, -2)] = e1*sigma_0
dn_lead[lat.neighbors()] = t0*sigma_0
syst.attach_lead(dn_lead)

prim_vecs=tinyarray.array([(a,0.,0.),(0.,a,0.),(0.,0.,a)])
offset=tinyarray.array((-1.2663118961, 0.4114496766,0.0))
lat2=kwant.lattice.Monatomic(prim_vecs, offset, norbs=no)

syst[lat2(0, 0, N)] = e1*sigma_0
syst[amorphous_lat(2*N-1), lat2(0, 0, N)] = tr*sigma_0

up_lead = kwant.Builder(kwant.TranslationalSymmetry((0, 0, a)))
up_lead[lat2(0, 0, N+1)] = e1*sigma_0
up_lead[lat2.neighbors()] = t0*sigma_0
syst.attach_lead(up_lead)

system=kwant.plot(syst, site_lw=0.1, site_color=family_color,
hop_lw=hopping_lw)


trans=True
if trans:
syst = syst.finalized()
energies = []
data = []

for ie in range(-320,520):
energy = ie * 0.001
smatrix = kwant.smatrix(syst, energy=energy)
energies.append(energy)
data.append(0.5*smatrix.transmission(1, 0))
fig = pyplot.figure(figsize=(6,2))
pyplot.plot(energies, data)
pyplot.xlim([-0.32,0.52])
pyplot.ylim([-0.03,1.25])
pyplot.xlabel("energy [eV]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()


strans1=True
if strans1:
#syst = syst.finalized()
energies = []
data = []

def oscle(ene, lead_nr):
wfs=kwant.wave_function(syst, ene,
check_hermiticity=True)(lead_nr)
spin_current_z = 0
for psi in wfs:
psi_start = psi[0 : 2]
psi_end = psi[2 * 61: 2 * 61 + 2]
#spin_current_z += -2 *
imag(psi_end.conjugate().dot(sigma_z).dot(psi_start))
spin_current_z +=
abs(imag(psi_end.conjugate().dot(sigma_z).dot(psi_start)))
return spin_current_z
for ie in range(-320,520):
energy = ie * 0.001
energies.append(energy)
data.append(oscle(ene=energy, lead_nr=1))
fig = pyplot.figure(figsize=(6,2))
pyplot.plot(energies, data)
pyplot.show()




strans=True
if strans:
#syst = syst.finalized()
J_spin = kwant.operator.Current(syst, sigma_z,
where=[(amorphous_lat(0), amorphous_lat(N-2))], sum=True)
all_wfs = kwant.wave_function(syst, energy=0.25)(1)
spin_current_list = sum(J_spin(wf) for wf in all_wfs)
print(spin_current_list)

[Joseph Weston]

Hi Patrik,

Post by Patrik Arvoy
1- The wave function approach does not give the the up/down
conductance similar to total.
psi_start = psi[0 : 2]
psi_end = psi[2 * 61: 2 * 61 + 2]

Do you know if the first and last sites sites as they are ordered in the
finalized system are indeed the sites that you want? In your script you
never refer to the 'sites' attribute of the finalized system, so I would
guess that these are not the sites that you want; sites number 0 and 61
are probably not even connected.

Post by Patrik Arvoy
2- the current-operator approach gives the error
'kwant.graph.core.EdgeDoesNotExistError'.

This means that you are trying to calculate the current between 2 sites
that are not connected by a hopping.

Post by Patrik Arvoy
J_spin = kwant.operator.Current(syst, sigma_z,
where=[(amorphous_lat(0), amorphous_lat(N-2))], sum=True)

I see that you are trying to calculate the current between
'amorphous_lat(0)' and 'amorphous_lat(N-2)', and if I colour these sites
when plotting the system I see indeed that they are at opposite ends of
your double helix structure, and that they are not connected by a hopping.


Happy Kwanting,


Joe

[Patrik Arvoy]

Hi Joe,

Thank you for the reply!

As you mentioned they are not connected via a single hopping. I was trying
to get up/down conductance through the system
from one end to the other, similar to the conductance from lead 0 to lead 1.

What would be the equivalent of conductance from lead 0 to lead 1 (in the
'if trans'), but for either up or down spin?

Regards
Patrik

Post by Joseph Weston
Hi Patrik,

Post by Patrik Arvoy
1- The wave function approach does not give the the up/down
conductance similar to total.
psi_start = psi[0 : 2]
psi_end = psi[2 * 61: 2 * 61 + 2]

Do you know if the first and last sites sites as they are ordered in the
finalized system are indeed the sites that you want? In your script you
never refer to the 'sites' attribute of the finalized system, so I would
guess that these are not the sites that you want; sites number 0 and 61
are probably not even connected.

Post by Patrik Arvoy
2- the current-operator approach gives the error
'kwant.graph.core.EdgeDoesNotExistError'.

This means that you are trying to calculate the current between 2 sites
that are not connected by a hopping.

Post by Patrik Arvoy
J_spin = kwant.operator.Current(syst, sigma_z,
where=[(amorphous_lat(0), amorphous_lat(N-2))], sum=True)

I see that you are trying to calculate the current between
'amorphous_lat(0)' and 'amorphous_lat(N-2)', and if I colour these sites
when plotting the system I see indeed that they are at opposite ends of
your double helix structure, and that they are not connected by a hopping.
Happy Kwanting,
Joe

[Joseph Weston]

Hi Patrik,

I was trying to get up/down conductance through the system
from one end to the other, similar to the conductance from lead 0 to lead 1.
What would be the equivalent of conductance from lead 0 to lead 1 (in
the 'if trans'), but for either up or down spin?

Ah, this is a new feature introduced in Kwant 1.3! The way to do this is
to first tell Kwant that your leads have a conservation law associated
with them, and to tell Kwant the basis that you want to use for the
modes in the lead.
Without this information the modes Kwant calculates in the leads will be
an arbitrary superposition of spin up and spin down, and you will not
easily be able to separate these.

You provide this information with the 'conservation_law' parameter when
creating the Builder for your leads. This parameter should be a
Hermitian matrix with integer eigenvalues, who's eigenvectors you want
to take as the basis for the modes. The modes will be ordered by their
eigenvalues. This is explained in the docs[1]; the example given is for
superconductivity, but the principle is the same. For your case a
reasonable choice would be sigma_z:

sym = kwant.TranslationalSymmetry([0, 0, -a])
dn_lead = kwant.Builder(sym, conservation_law=-sigma_z)



There is a minus sign in front of the sigma_z so that the "up" modes
come first, and then the "down" modes.

Then, when computing the transmission, instead of saying

smatrix.transmission(1, 0)

which will give you the total transmission from lead 0 to lead 1, you
can say

smatrix.transmission((1, 0), (0, 0))

to get the transmission from the spin up block (block 0) of lead 0 to
the spin up block of lead 1. The first number in each pair is the lead
number, and the second is the block index with respect to the
conservation law you defined when creating the Builder (0 being spin up,
and 1 being spin down in this case).

Similarly you can use

smatrix.transmission((1, 1), (0, 1))

for transmission from spin down to spin down, or any other combination
to calculate the transmissions between different spin blocks.

Hope that helps,

Joe

[1]: https://kwant-project.org/doc/1/tutorial/superconductors

[Patrik Arvoy]

Hi Joe,

Perfect, it works now.
Thanks a lot!

Regards
Patrik

Post by Joseph Weston
Hi Patrik,

I was trying to get up/down conductance through the system
from one end to the other, similar to the conductance from lead 0 to lead 1.
What would be the equivalent of conductance from lead 0 to lead 1 (in
the 'if trans'), but for either up or down spin?

Ah, this is a new feature introduced in Kwant 1.3! The way to do this is
to first tell Kwant that your leads have a conservation law associated
with them, and to tell Kwant the basis that you want to use for the
modes in the lead.
Without this information the modes Kwant calculates in the leads will be
an arbitrary superposition of spin up and spin down, and you will not
easily be able to separate these.
You provide this information with the 'conservation_law' parameter when
creating the Builder for your leads. This parameter should be a
Hermitian matrix with integer eigenvalues, who's eigenvectors you want
to take as the basis for the modes. The modes will be ordered by their
eigenvalues. This is explained in the docs[1]; the example given is for
superconductivity, but the principle is the same. For your case a
sym = kwant.TranslationalSymmetry([0, 0, -a])
dn_lead = kwant.Builder(sym, conservation_law=-sigma_z)


There is a minus sign in front of the sigma_z so that the "up" modes
come first, and then the "down" modes.
Then, when computing the transmission, instead of saying
smatrix.transmission(1, 0)
which will give you the total transmission from lead 0 to lead 1, you
can say
smatrix.transmission((1, 0), (0, 0))
to get the transmission from the spin up block (block 0) of lead 0 to
the spin up block of lead 1. The first number in each pair is the lead
number, and the second is the block index with respect to the
conservation law you defined when creating the Builder (0 being spin up,
and 1 being spin down in this case).
Similarly you can use
smatrix.transmission((1, 1), (0, 1))
for transmission from spin down to spin down, or any other combination
to calculate the transmissions between different spin blocks.
Hope that helps,
Joe
[1]: https://kwant-project.org/doc/1/tutorial/superconductors