-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathLandauer_v56_B_15.m
191 lines (160 loc) · 6.96 KB
/
Landauer_v56_B_15.m
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
%% Version 56-B-15: This code is to calculate dI/dV false plot/peak spacings as a function of ng and Vz for the multi-level system.
% Soft gap is included in the topological regime.
clear;
tic;
C = 20; part = 15;
%% Parameters Setting
% Note that the length scale is in unit of lattice constant, which is 10nm.
% Basic Parameters
N = 20; % The minimum of N should be chosen as the # of energy levels we consider. i.e. N >= levelN.
vec_N = [(N-1) N (N+1)];
t = 25; %unit: meV
Delta_0 = 0.9; %unit: meV
Vzc = 4.2;
wireLength = 150; %unit: 10nm
alpha = 2.5; %unit: meV
mu = 2.5; %unit: meV
lambda = 1.4; %unit: meV
T = 0.01; %unit: meV
Ec = 3; %unit: meV
VD1 = 1; %unit: meV
VD2 = 4; %unit: meV
N_dot = 26; %unit: 10nm
Nbarrier = 2; %unit: 10nm
Ebarrier = 10; %unit: meV
N_tot = wireLength;
s = 1e-3; % Vstep = 1e-3; resolution is 5 times smaller than the width.
rho_F = 0.1;
V_SC = 10^5;
V_metal = 2.5;
threshold = 1e-2;
%% Construct the variable space
VzMin = 0; VzMax = 5.5; VzNumber = 551;
VzStep = (VzMax - VzMin)./(VzNumber - 1);
VzRange = linspace(VzMin,VzMax,VzNumber);
Vc2 = sqrt(lambda.^2 + mu.^2);
Vc2Number = floor(1 + (Vc2 - VzMin)./VzStep);
ngMin = 19; ngMax = 21; ngNumber = 2001; % (19.4; 19.6; 81)
ngStep = (ngMax - ngMin)./(ngNumber - 1);
ngRange = linspace(ngMin,ngMax,ngNumber);
Vmin = 0; Vmax = 1.1; Vnumber = 11001; % (-0.0025, 4, 1002)1101
Vstep = (Vmax - Vmin)./(Vnumber - 1); % = 1e-3
Vrange = linspace(Vmin,Vmax,Vnumber);
G_even = zeros(VzNumber,ngNumber);
G_odd = zeros(VzNumber,ngNumber);
G_odd_L = zeros(VzNumber,ngNumber);
G_odd_R = zeros(VzNumber,ngNumber);
G_even_L = zeros(VzNumber,ngNumber);
G_even_R = zeros(VzNumber,ngNumber);
level_num = zeros(1,C);
total_num = zeros(1,C);
E_less_T_num = zeros(1,C);
Q_0 = zeros(1,C);
Parity = zeros(1,C);
%%
parfor k = 1:C % VzNumber
K = C.*part+k;
%K = C.*part+k-10;
Vz = VzMin + (K-1).*VzStep;
disp(K);
Delta1 = Delta_0.*sqrt(1 - (Vz./Vzc).^2).*(Vz<Vzc);
%Delta1 = Delta_0;
Delta2 = Delta_0.*sqrt(1 - (Vz./Vc2).^2).*(Vz<Vc2);
Gamma_L = [];
Gamma_R = [];
Lambda_L = [];
Lambda_R = [];
vecE_temp = [];
parity = 0;
Q_0(k) = 1;
tempE = Delta2 + 1e-16;
metal2_E = [tempE];
accuD = 0;
n1 = 0;
while tempE < V_metal
n1 = n1+1;
rho_w = dosH_se_v16(t,Delta1,Delta2,N_tot,alpha,mu,VD1,VD2,N_dot,Nbarrier,Ebarrier,Vz,lambda,tempE,s);
[dege_WF,lambda_m] = eigs(rho_w,4.*N_tot,'LM');
lambda_m = diag(lambda_m);
lambda_m = real(lambda_m);
index = find(lambda_m >= threshold);
d = length(index);
accuD = accuD + d;
F = accuD./V_SC;
if (tempE <= Delta1) && (tempE > Delta2)
dEn = sqrt(Delta2.^2 + (accuD./(2.*rho_F.*V_SC)).^2) - tempE;
elseif (tempE > Delta1)
dEn = sqrt((Delta1.^2 + Delta2.^2 + (F./rho_F).^2).^2 - 4.*(Delta1.*Delta2).^2)./(2.*F./rho_F) - tempE;
end
lambda_m = lambda_m(index); % selected d degeneracy
dege_WF = dege_WF(:,index); % selected d degeneracy
prefactor = sqrt(lambda_m'.*dEn); % (4.*N_tot)-by-1 vecotr
dege_WF = prefactor.*dege_WF;
vecE_temp = [vecE_temp repmat(tempE,1,d)];
tempE = tempE + dEn;
metal2_E = [metal2_E tempE];
Gamma2_l_temp = not(parity).*(abs(dege_WF(4.*Nbarrier+3,:)).^2 + abs(dege_WF(4.*Nbarrier+4,:)).^2) + parity.*(abs(dege_WF(4.*Nbarrier+1,:)).^2 + abs(dege_WF(4.*Nbarrier+2,:)).^2);
Gamma2_r_temp = not(parity).*(abs(dege_WF(4*N_tot - 1,:)).^2 + abs(dege_WF(4*N_tot,:)).^2) + parity.*(abs(dege_WF(4*N_tot - 3,:)).^2 + abs(dege_WF(4*N_tot - 2,:)).^2);
Lambda2_l_temp = not(parity).*(abs(dege_WF(4.*Nbarrier+1,:)).^2 + abs(dege_WF(4.*Nbarrier+2,:)).^2) + parity.*(abs(dege_WF(4.*Nbarrier+3,:)).^2 + abs(dege_WF(4.*Nbarrier+4,:)).^2);
Lambda2_r_temp = not(parity).*(abs(dege_WF(4*N_tot - 3,:)).^2 + abs(dege_WF(4*N_tot - 2,:)).^2) + parity.*(abs(dege_WF(4*N_tot - 1,:)).^2 + abs(dege_WF(4*N_tot,:)).^2);
Gamma_L = [Gamma_L Gamma2_l_temp];
Gamma_R = [Gamma_R Gamma2_r_temp];
Lambda_L = [Lambda_L Lambda2_l_temp];
Lambda_R = [Lambda_R Lambda2_r_temp];
end
level_num(k) = length(metal2_E)-1;
vecE = vecE_temp;
temp = tanh(vecE./(2.*T));
total_num(k) = length(vecE);
Total_N = total_num(k);
y = zeros(1,Total_N);
x1 = vecE(1)./(2.*T);
y(1) = 2.*x1 + log((1+exp(-2.*x1))./2);
for p = 2:Total_N
x = vecE(p)./(2.*T);
z = 2.*x + log((1+exp(-2.*x))./2);
a = max([y(p-1) z]);
b = min([y(p-1) z]);
y(p) = b - log(1 + exp(-(a-b)) - exp(-a));
end
gamma = exp(-y(Total_N));
prod_T = 1-gamma;
prod_T_noP = prod_T./temp;
temp_index = find(vecE<T);
E_less_T_num(k) = length(temp_index);
for m = 1:ngNumber
ng = ngMin + (m-1).*ngStep;
U_N = Ec.*(vec_N-ng).^2;
dU = diff(U_N); % dU = U(N) - U(N-1)
if mod(N,2)==0 % Make dU(1) is for N odd and dU(2) is for N even.
dU = fliplr(dU);
end
Z_odd = (1 + exp(-vecE./T)).*(gamma + exp(-(-Q_0(k).*dU(1)./T)).*(2-gamma));
Z_even = (1 + exp(-vecE./T)).*(gamma + exp(-(Q_0(k).*dU(2)./T)).*(2-gamma));
tempM_odd_L = zeros(1,total_num(k));
tempM_odd_R = zeros(1,total_num(k));
tempM_even_L = zeros(1,total_num(k));
tempM_even_R = zeros(1,total_num(k));
for n = 0:1
for Q = -1:2:1 % Q is either (-1) or (+1)
dE_odd = vecE.*(1 - 2.*n) + Q.*dU(1); % N is odd.
dE_even = vecE.*(1 - 2.*n) - Q.*dU(2); % N is even.
Fp_odd = 0.5.*exp(-(vecE - Q_0(k).*dU(1))./(2.*T)).*sech(dE_odd./(2.*T))./T.*(1 + Q.*Q_0(k).*((-1).^n).*prod_T_noP)./Z_odd;
Fp_even = 0.5.*exp(-(vecE + Q_0(k).*dU(2))./(2.*T)).*sech(dE_even./(2.*T))./T.*(1 + Q.*Q_0(k).*((-1).^n).*prod_T_noP)./Z_even;
tempM_odd_L = tempM_odd_L + Fp_odd.*((1-n).*Gamma_L + n.*Lambda_L);
tempM_odd_R = tempM_odd_R + Fp_odd.*((1-n).*Gamma_R + n.*Lambda_R);
tempM_even_L = tempM_even_L + Fp_even.*((1-n).*Gamma_L + n.*Lambda_L);
tempM_even_R = tempM_even_R + Fp_even.*((1-n).*Gamma_R + n.*Lambda_R);
end
end
G_odd_L(k,m) = sum(tempM_odd_L); % Sum over the changed energy states and configurations
G_odd_R(k,m) = sum(tempM_odd_R); % Sum over the changed energy states and configurations
G_even_L(k,m) = sum(tempM_even_L); % Sum over the changed energy states and configurations
G_even_R(k,m) = sum(tempM_even_R); % Sum over the changed energy states and configurations
G_odd(k,m) = G_odd_L(k,m).*G_odd_R(k,m)./(G_odd_L(k,m) + G_odd_R(k,m));
G_even(k,m) = G_even_L(k,m).*G_even_R(k,m)./(G_even_L(k,m) + G_even_R(k,m));
disp(['(Vz,ng)=(',num2str(Vz),',',num2str(ng),')']);
end
end
toc;
save(['Landauer_v56_L=150_partB_',num2str(part),'.mat'])