-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmezo.py
205 lines (166 loc) · 6.86 KB
/
mezo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
import numpy
import matplotlib
from matplotlib import pylab, mlab, pyplot
from pylab import *
from numpy import *
from numpy.linalg import *
import scipy.linalg as sl # we need a generailzed eigen solver
class lead1D:
'A class for simple 1D leads'
def __init__(self,eps0=0,gamma=-1,**kwargs):
'We assume real hopping \gamma and onsite \epsilon_0 parameters!'
self.eps0=eps0
self.gamma=gamma
return
def Ek(self,k):
'Spectrum as a function of k'
return self.eps0+2*self.gamma*cos(k)
def kE(self,E,**kwargs):
'''
Spectrum as a function of E.
If keyword a=True is given than
it gives back two k values,
one positive and one negative.
'''
a = kwargs.get('a',False)
k=arccos((E-self.eps0)/(2*self.gamma))
if a:
return array([-k,k])
else:
return k
def vE(self,E=0,**kwargs):
'''
Group velocity as a function of E.
If keyword a=True is given than
it gives back two v values,
one positive and one negative.
'''
a = kwargs.get('a',False)
k=self.kE(E)
v= -2*self.gamma*sin(k)
if a:
return array([-v,v])
else:
return v
def sgf(self,E=0):
'''
Surgace Green's function of a seminfinite 1D lead.
'''
return exp(1.0j *self.kE(E))/self.gamma
def sgfk(self,k=pi/2):
'''
Surgace Green's function of a seminfinite 1D lead in terms of k.
'''
return exp(1.0j*k)/self.gamma
def vk(self,k=pi/2):
'''
Group velocity in terms of k
'''
return -2*self.gamma*sin(k)
class lead:
'''
This is a class for a generic lead.
Without degenerate subspace ironing!
Without
'''
def __init__(self,H0,H1,**kwargs):
'''
Initialization of the leads.
Maybe some checks could come here later on
or even the svd...
'''
self.H0=matrix(H0)
self.H1=matrix(H1)
self.dim=shape(H0)[0]
self.tol=1.0e-10
def get_spectrEk(self,kran,**kwargs):
'''Extract spectrum as function of momentum E(k)'''
spectr=zeros((len(kran),self.dim))
for i in range(len(kran)):
k=kran[i]
spectr[i,:]=eigvalsh(self.H0+exp(1.0j*k)*self.H1+exp(-1.0j*k)*self.H1.H)
return spectr
def set_ene(self,E,**kwargs):
'''
Build energy dependent quantities
This is the part where the bulk and surface
Green's function and self-energies are calculated.
'''
self.E=E # This is the enregy we are using
Z=zeros_like(self.H0) # Zero matrix and Identity matrix
Id=eye(self.dim) # of size dimXdim
# Setting up the tricky eigenvalue problem ...
A =vstack((hstack((E*eye(self.dim)-self.H0,-self.H1.H)),
hstack((Id,Z))))
B =vstack((hstack((self.H1,Z)),
hstack((Z,Id))))
# and solving it.
w,vec=sl.eig(A,B)
self.w=w
# get k if asked for
if kwargs.get('get_k'):
self.k=-1.0j*log(w)
# Obtaining and normailizing the eigen vectors of the transverse modes
vec=matrix(vec[:self.dim,:])
vg=(zeros_like(w));
for i in range(2*self.dim):
vec[:,i]=vec[:,i]/norm(vec[:,i])
vg[i]=1.0j*((vec[:,i].H*(w[i]*self.H1-w[i]**(-1)*self.H1.H)*vec[:,i])[0,0])
self.vg=vg
# Sorting left and right channels in to propagaintg and decaying
opened=where((abs(abs(w)-1)<self.tol))
closed=where(abs(abs(w)-1)>self.tol)
left =union1d(intersect1d(opened,where(real(vg)<0)),
intersect1d(closed,where(abs(w)>1)))
right=union1d(intersect1d(opened,where(real(vg)>0)),
intersect1d(closed,where(abs(w)<1)))
self.vec_left=vec[:,left]
self.w_left =w[left]
self.vg_left =vg[left]
self.vec_right=vec[:,right]
self.w_right =w[right]
self.vg_right =vg[right]
if len(self.w_right)!=len(self.w_left):
print('Problem with partitioning!!')
return
# Obtaining the duals
self.vec_right_dual=inv(self.vec_right)
self.vec_left_dual=inv(self.vec_left)
# Packaging open channel quantities
left_open =where(abs(abs(w[left])-1 )<self.tol)
right_open =where(abs(abs(w[right])-1 )<self.tol)
self.vg_left_open=(self.vg[left])[left_open]
self.vg_right_open=(self.vg[right])[right_open]
self.w_left_open=(self.w[left])[left_open]
self.w_right_open=(self.w[right])[right_open]
self.vec_left_open =(self.vec_left[:,left_open])[:,0,:] # for some wierd reson the [...] is needed
self.vec_right_open =(self.vec_right[:,right_open])[:,0,:] # if i want to have a nice matrix
self.vec_left_dual_open =(self.vec_left_dual[left_open,:])[0,:,:] # same here with the [...]
self.vec_right_dual_open=(self.vec_right_dual[right_open,:])[0,:,:]
# Get transfere matrices
self.T_leftm1 = self.vec_left * matrix(diag(w[left])) * self.vec_left_dual
self.T_left = self.vec_left * matrix(diag(1/w[left])) * self.vec_left_dual
self.T_right = self.vec_right * matrix(diag(w[right])) * self.vec_right_dual
self.T_rightm1 = self.vec_right * matrix(diag(1/w[right])) * self.vec_right_dual
self.V=self.H1*(self.T_leftm1-self.T_right)
if kwargs.get('get_selfE'):
# Get self-energy
self.selfEL=self.H1.H*self.T_left
self.selfER=self.H1*self.T_right
elif kwargs.get('get_bulkG'):
# Get Green's functions of infinite bulk and leads
self.g00=inv(self.V)
self.gsL=self.T_left*inv(self.H1)
self.gsR=self.T_right*inv(self.H1.H)
else:
self.gsL=self.T_left*inv(self.H1)
self.gsR=self.T_right*inv(self.H1.H)
def Tz(self,z,**kwargs):
'''
A function to obtain arbitrary tranfere matrices
'''
if (hasattr(self,'E')):
return (self.vec_left * matrix(diag(self.w_left**z)) * self.vec_left_dual,
self.vec_right * matrix(diag(self.w_right**z)) * self.vec_right_dual)
else:
print('Set energy first')