Dear Sir, I was trying to compute the transmissions for a 2 terminal GNR from lead 0 (up/down) to lead 1 (up/down) using the following discussion. http://kwant-discuss.kwant-project.narkive.com/ZRoAKlBg/up-down-conductance I constructed the leads as follows with the help of tutorial of superconductor: def create_lead_h(W, symmetry, axis=(0, 0)): lead = kwant.Builder(symmetry,conservation_law=-sigma_z) lead[lat.wire(axis, W)] = 0. * sigma_z lead[lat.neighbors(1)] = -1.*sigma_z return lead and have calculated the following transmissions: smatrix = kwant.smatrix(sys, e, args=(params,)) su = smatrix.transmission((1, 0), (0, 0)) + smatrix.transmission((1, 0), (0, 1)) sd = smatrix.transmission((1, 1), (0, 0)) + smatrix.transmission((1, 1), (0, 1)) cc = smatrix.transmission(1, 0) The system has Rashba interaction. My concern is when I turn off the Rashba interaction, su and sd should be equal to each other and have integer values. But I'm getting the integer values only for su. I don't understand why I'm getting this behaviour. The program is written here for better understanding. ================================================================= from math import sqrt import random import itertools as it import tinyarray as ta import numpy as np from matplotlib import pyplot import kwant #from sympy import S, Eq, solve, Symbol class Honeycomb(kwant.lattice.Polyatomic): """Honeycomb lattice with methods for dealing with hexagons""" def __init__(self, name=''): prim_vecs = [[0.5, sqrt(3)/2], [1, 0]] # bravais lattice vectors # offset the lattice so that it is symmetric around x and y axes basis_vecs = [[-0.5, -1/sqrt(12)], [-0.5, 1/sqrt(12)]] super(Honeycomb, self).__init__(prim_vecs, basis_vecs, name, norbs=2) self.a, self.b = self.sublattices def hexagon(self, tag): """ Get sites belonging to hexagon with the given tag. Returns sites in counter-clockwise order starting from the lower-left site. """ tag = ta.array(tag) # a-sites b-sites deltas = [(0, 0), (0, 0), (1, 0), (0, 1), (0, 1), (-1, 1)] lats = it.cycle(self.sublattices) return (lat(*(tag + delta)) for lat, delta in zip(lats, deltas)) def hexagon_neighbors(self, tag, n=1): """ Get n'th nearest neighbor hoppings within the hexagon with the given tag. """ hex_sites = list(self.hexagon(tag)) return ((hex_sites[(i+n)%6], hex_sites[i%6]) for i in range(6)) def cross(W, L): def shape(pos): return ((-W <= pos[1] <= W and -L <= pos[0] <= L)) # horizontal return shape def circle(W,L,r1,r2): def ring(pos): (x, y) = pos rsq = x ** 2 + y ** 2 if (sqrt(rsq)>=r1): return (-W <= pos[1] <= W and -L <= pos[0] <= L) return ring ## Pauli matrices ## sigma_0 = np.array([[1, 0], [0, 1]]) # identity sigma_x = np.array([[0, 1], [1, 0]]) sigma_y = np.array([[0, -1j], [1j, 0]]) sigma_z = np.array([[1, 0], [0, -1]]) ## Hamiltonian value functions ## def onsite_potential(site, params): return params['ep'] * sigma_0 def potential_shift(site, params): return -params['mu'] * sigma_0 def kinetic(site_i, site_j, params): return -params['gamma'] * sigma_0 def rashba(site_i, site_j, params): d_ij = site_j.pos - site_i.pos rashba = 1j * params['V_R'] * (sigma_x * d_ij[1] - sigma_y * d_ij[0]) return rashba + kinetic(site_i, site_j, params) def spin_orbit(site_i, site_j, params): so = 1j * params['V_I'] * sigma_z return so lat = Honeycomb() pv1, pv2 = lat.prim_vecs ysym1 = kwant.TranslationalSymmetry(pv2 - 2*pv1) # lattice symmetry in -y direction ysym2 = kwant.TranslationalSymmetry(-pv2 + 2*pv1) # lattice symmetry in -y direction xsym1 = kwant.TranslationalSymmetry((-pv2)) # lattice symmetry in -x direction xsym1.add_site_family(lat.sublattices[0], other_vectors=[(-2, 1)]) xsym1.add_site_family(lat.sublattices[1], other_vectors=[(-2, 1)]) xsym2=kwant.TranslationalSymmetry((pv2)) xsym2.add_site_family(lat.sublattices[0], other_vectors=[(-2, 1)]) xsym2.add_site_family(lat.sublattices[1], other_vectors=[(-2, 1)]) def create_lead_h(W, symmetry, axis=(0, 0)): lead = kwant.Builder(symmetry,conservation_law=-sigma_z) lead[lat.wire(axis, W)] = 0. * sigma_z lead[lat.neighbors(1)] = -1.*sigma_z#kinetic return lead def create_lead_v(L,W, symmetry): axis=(0, 0) lead = kwant.Builder(symmetry,conservation_law=-sigma_z) lead[lat.wire(axis, W)] = 0. * sigma_0 lead[lat.neighbors(1)] = kinetic return lead def create_hall_cross(W, L): ## scattering region ## sys = kwant.Builder() sys[lat.shape(circle(W,L,r1,r2), (0, r1+1))] = onsite_potential sys[lat.neighbors(1)] = rashba sys[lat.neighbors(2)] = spin_orbit ## leads ## leads = [create_lead_h(W, xsym1)] leads += [create_lead_h(W, xsym2 for lead in leads: sys.attach_lead(lead) return sys #========================================================================= def spin_conductance(G, lead_out, lead_in, sigma): """Calculate the spin conductance between two leads. Parameters ---------- G : an instance of `kwant.solvers.common.GreensFunction` The Greens function of the system as returned by `kwant.greens_function`. lead_out : integer The lead where spin current is collected lead_in : integer The lead where spin current is injected sigma : `numpy.ndarray` of shape (2, 2) The Pauli matrix of the quantization axis along which to measure the spin current Notes ----- Calculates the spin conductance between two leads p and q according to: G_{pq} = Tr[ σ_{α} Γ_{q} G_{qp} Γ_{p} G^+_{qp} ] Where Γ_{q} is the coupling matrix to lead q ( = i[Σ - Σ^+] ) and G_{qp} is the submatrix of the system's Greens function between sites interfacing with lead p and q. """ ttdagger = G._a_ttdagger_a_inv(lead_out, lead_in) sigma_z_matrix = np.kron(np.eye(ttdagger.shape[0]//2), sigma) return np.trace(sigma_z_matrix.dot(ttdagger)).real#, np.trace(ttdagger).real if __name__ == '__main__': params = dict(gamma=1., ep=0, mu=0.0, V_I=0.0, V_R=0.0) W=17.5 L=20 r1=0 r2=10 sys = create_hall_cross(W, L) kwant.plot(sys, site_color=site_color, site_size=site_size) sys=sys.finalized() band=100 Es = np.linspace(-0.5, 0.5, band) sp_up=[] sp_down=[] ccurrent=[] spz=[] for e in Es: G = kwant.greens_function(sys, e, args=(params,)) smatrix = kwant.smatrix(sys, e, args=(params,)) su = smatrix.transmission((1, 0), (0, 0)) + smatrix.transmission((1, 0), (0, 1)) sd = smatrix.transmission((1, 1), (0, 0)) + smatrix.transmission((1, 1), (0, 1)) cc = smatrix.transmission(1, 0) spin_pol_z = spin_conductance(G, 1, 0, sigma_z) diff=su-sd print (e,su,sd,cc,spin_pol_z,diff) sp_up.append(su) sp_down.append(sd) ccurrent.append(cc) spz.append(spin_pol_z) pyplot.figure() pyplot.plot(Es, sp_up) pyplot.plot(Es, sp_down) pyplot.plot(Es, ccurrent) pyplot.xlabel("energy [t]") pyplot.ylabel("spin polarization_x [e^2/h]") pyplot.show() ============================================================= With Regards, Sudin