import kwant
import numpy as np
class WavefunctionLead(kwant.builder.InfiniteSystem):
def __init__(self, lead):
self.__dict__ = lead.__dict__
def modes(self, energy, args=()):
prop_modes_fix, stab_modes_fix = \
super(WavefunctionLead, self).modes(energy=0.4, args=args)
size_fix = prop_modes_fix.wave_functions.shape[1]
fix = []
for i in range(size_fix):
fix0 = prop_modes_fix.wave_functions[0][i]
fix.append(fix0)
prop_modes, stab_modes = \
super(WavefunctionLead, self).modes(energy=energy, args=args)
wf_size = prop_modes.wave_functions.shape[0] #number of orbital
num_modes = prop_modes.wave_functions.shape[1]#number of modes (2*num_modes)
## change of the phase
for i in range(wf_size):
for j in range(num_modes):
f = fix[j]
prop_modes.wave_functions[i,j]=prop_modes.wave_functions[i,j]*(prop_modes.wave_functions[i,j]/f)
if prop_modes.wave_functions[i,j].real >= 0 :
prop_modes.wave_functions[i,j] = -1 * prop_modes.wave_functions[i,j]
else :
continue
return prop_modes, stab_modes
def make_system(a, t, L, f, r, rcc, kso, w):
lat = kwant.lattice.square(a)
sys = kwant.Builder()
def myShape(pos):
(x, y) = pos
if (L/2-f<= y < L/2+f) and (0<=x <2):
return True
elif (L/2-f<= y < L/2+f) and (L-1<= x< L+1):
return True
elif ( 2 <= x < L-1 ) and ( 0 <= y < L-1) and rcc ** 2 <(x ** 2 + y ** 2) and (3*r ** 2 < ((x-(L/2)) ** 2 + (y-(L/2)) ** 2)):
return True
def onsite(site1, tm):## Electric field in x direction
x1 =site1.pos[0]
y1 = site1.pos[1]
return 4*t + np.sin(w*tm)*w*x1*kso
sys[lat.shape(myShape, (0, L/2))] = onsite
sys[kwant.builder.HoppingKind((1, 0), lat, lat)] = -t
sys[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t
##leads##
lead0 = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
def lead_shape_0(pos):
(x, y) = pos
return ( L/2-f <= y < L/2+f)
lead0[lat.shape(lead_shape_0, (0, L/2))] = 4*t
lead0[kwant.builder.HoppingKind((1, 0), lat, lat)] = -t
lead0[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t
lead1 = kwant.Builder(kwant.TranslationalSymmetry((a, 0)))
def lead_shape_1(pos):
(x, y) = pos
return ( L/2-f <= y < L/2+f)
lead1[lat.shape(lead_shape_1, (0, L/2))] = 4*t
lead1[kwant.builder.HoppingKind((1, 0), lat, lat)] = -t
lead1[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t
sys.attach_lead(lead0)
sys.attach_lead(lead1)
fsys = sys.finalized()
###
fsys.leads = [WavefunctionLead(lead) for lead in fsys.leads]
###
return fsys
def main():
## w = frequency
## kso = Rashba spin orbit coupling
sys = make_system(a=1, t=1, L=50, f=2, r=6, rcc=5, kso = 0.3/50, w = 0.001)
energies=[0.00001 * i+0.40 for i in xrange(300)]
tm =1000 #time parameter
for energy in energies:
S = kwant.smatrix(sys, energy, args=[tm])
mode = S.lead_info
wf =mode[0].wave_functions
print wf
if __name__ == '__main__':
main()