Hi Nafise

You can access the s-matrix elements using the submatrix method (https://kwant-project.org/doc/latest/reference/generated/kwant.solvers.common.SMatrix).

Best,

Antonio

On Sun, 18 Apr 2021, 11:55 Nafise Nouri, <nafise.nour@gmail.com> wrote:

Dear all,
I need to access the matrix diagonal element in the scattering region in my program. Is there any way to access  the matrix diagonal element in the following program or do I have to use numpy?

Best regard,

import kwant
import matplotlib.pyplot as plt
import tinyarray
import numpy as np
import math
from functools import partial
from kwant.physics import dispersion

#mport matplotlib

d=1.42;
a1=d*math.sqrt(3);

latt = kwant.lattice.general([(a1,0),(0.5*a1,a1*math.sqrt(3)/2)],
                             [(0,0),(0,a1/math.sqrt(3))])
a,b = latt.sublattices
syst= kwant.Builder()

#...................................................................................
def rectangle(pos):
    x, y = pos
    z=x**2+y**2
    return -4*3*a1<y<4.5*3*a1 and -4.5*d<=x<4.5*d

syst[a.shape(rectangle, (2,2))] = 0
syst[b.shape(rectangle, (2,2))] = 0

#nearest neighbors.............................................................
syst[[kwant.builder.HoppingKind((0,0),a,b)]] = -2.7
syst[[kwant.builder.HoppingKind((0,1),a,b)]] = -2.7
syst[[kwant.builder.HoppingKind((-1,1),a,b)]] = -2.7

ax=kwant.plot(syst);

sym = kwant.TranslationalSymmetry(latt.vec((-1,0)))

sym.add_site_family(latt.sublattices[0], other_vectors=[(-1, 2)])
sym.add_site_family(latt.sublattices[1], other_vectors=[(-1, 2)])
lead = kwant.Builder(sym)


def lead_shape(pos):
    x, y = pos
    return  -4*3*a1<y<4.5*3*a1

lead[a.shape(rectangle, (2,2))] = 0
lead[b.shape(rectangle, (2,2))] = 0


#nearest neighbors.............................................................
lead[[kwant.builder.HoppingKind((0,0),a,b)]] = -2.7
lead[[kwant.builder.HoppingKind((0,1),a,b)]] = -2.7
lead[[kwant.builder.HoppingKind((-1,1),a,b)]] = -2.7

syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
ax=kwant.plot(syst);

###############################################################################\
fsys = syst.finalized()
Sites= list(fsys.sites)   #list of all the sites in the scattering region
number_of_sites=len(Sites)


bands = dispersion.Bands(fsys.leads[0])
momenta = np.linspace(-np.pi,np.pi, 100)
energies=[bands(k) for k in momenta]

np.savetxt('band1.txt',bands(-np.pi))


plt.plot(momenta,energies)

plt.xlabel(r'$K$')
plt.ylabel(r'$band structure (eV)$')
plt.ylim((-0.5,0.5))
plt.savefig('bandsarmBiStrain.pdf')
plt.show()