-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathSchroed3D_FEM_f.m
195 lines (172 loc) · 16.1 KB
/
Schroed3D_FEM_f.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
192
193
194
195
function[E,psi]=Schroed3D_FEM_f(x,y,z,V0,Mass,n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=6.62606896E-34; %% Planck constant [J.s]
hbar=h/(2*pi);
e=1.602176487E-19; %% electron charge [C]
m0=9.10938188E-31; %% electron mass [kg]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nx=length(x);
Ny=length(y);
Nz=length(z);
dx=x(2)-x(1);
dy=y(2)-y(1);
dz=z(2)-z(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Building of the operators %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Second derivative X %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DX2(Ny=3,Nx=3,Nz=3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% -2 0 0 | 1 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 -2 0 | 0 1 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 -2 | 0 0 1 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 1 0 0 |-2 0 0 | 1 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 1 0 | 0 -2 0 | 0 1 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 1 | 0 0 -2 | 0 0 1 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 1 0 0 |-2 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 1 0 | 0 -2 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 1 | 0 0 -2 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% =========================================================================================================== %
% 0 0 0 | 0 0 0 | 0 0 0 ||-2 0 0 | 1 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 -2 0 | 0 1 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 -2 | 0 0 1 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 1 0 0 |-2 0 0 | 1 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 1 0 | 0 -2 0 | 0 1 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 1 | 0 0 -2 | 0 0 1 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 1 0 0 |-2 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 1 0 | 0 -2 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 1 | 0 0 -2 || 0 0 0 | 0 0 0 | 0 0 0 %
% =========================================================================================================== %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 ||-2 0 0 | 1 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 -2 0 | 0 1 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 -2 | 0 0 1 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 1 0 0 |-2 0 0 | 1 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 1 0 | 0 -2 0 | 0 1 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 1 | 0 0 -2 | 0 0 1 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 1 0 0 |-2 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 1 0 | 0 -2 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 1 | 0 0 -2 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Axy = [ ones(1,(Nx-1)*Ny) zeros(1,Ny) ];
Axyz = [ repmat(Axy,1,Nz-1) ones(1,(Nx-1)*Ny) ];
DX2 = (-2)*diag(ones(1,Nx*Ny*Nz)) + (1)*diag(Axyz,-Ny) + (1)*diag(Axyz,Ny);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Second derivative Y %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DY2(Ny=3,Nx=3,Nz=3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% -2 1 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 1 -2 1 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 1 -2 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 |-2 1 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 1 -2 1 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 1 -2 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 |-2 1 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 1 -2 1 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 1 -2 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% =========================================================================================================== %
% 0 0 0 | 0 0 0 | 0 0 0 ||-2 1 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 1 -2 1 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 1 -2 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 |-2 1 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 1 -2 1 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 1 -2 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 |-2 1 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 1 -2 1 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 1 -2 || 0 0 0 | 0 0 0 | 0 0 0 %
% =========================================================================================================== %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 ||-2 1 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 1 -2 1 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 1 -2 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 |-2 1 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 1 -2 1 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 1 -2 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 |-2 1 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 1 -2 1 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 1 -2 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
BB=ones(1,Nx*Ny*Nz-1);
BB(Ny:Ny:end)=0;
DY2=(-2)*diag(ones(1,Nx*Ny*Nz)) + (1)*diag(BB,-1) + (1)*diag(BB,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Second derivative Z %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DZ2(Ny=3,Nx=3,Nz=3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% -2 0 0 | 0 0 0 | 0 0 0 || 1 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 -2 0 | 0 0 0 | 0 0 0 || 0 1 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 -2 | 0 0 0 | 0 0 0 || 0 0 1 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 |-2 0 0 | 0 0 0 || 0 0 0 | 1 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 -2 0 | 0 0 0 || 0 0 0 | 0 1 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 -2 | 0 0 0 || 0 0 0 | 0 0 1 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 |-2 0 0 || 0 0 0 | 0 0 0 | 1 0 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 -2 0 || 0 0 0 | 0 0 0 | 0 1 0 || 0 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 -2 || 0 0 0 | 0 0 0 | 0 0 1 || 0 0 0 | 0 0 0 | 0 0 0 %
% =========================================================================================================== %
% 1 0 0 | 0 0 0 | 0 0 0 ||-2 0 0 | 0 0 0 | 0 0 0 || 1 0 0 | 0 0 0 | 0 0 0 %
% 0 1 0 | 0 0 0 | 0 0 0 || 0 -2 0 | 0 0 0 | 0 0 0 || 0 1 0 | 0 0 0 | 0 0 0 %
% 0 0 1 | 0 0 0 | 0 0 0 || 0 0 -2 | 0 0 0 | 0 0 0 || 0 0 1 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 1 0 0 | 0 0 0 || 0 0 0 |-2 0 0 | 0 0 0 || 0 0 0 | 1 0 0 | 0 0 0 %
% 0 0 0 | 0 1 0 | 0 0 0 || 0 0 0 | 0 -2 0 | 0 0 0 || 0 0 0 | 0 1 0 | 0 0 0 %
% 0 0 0 | 0 0 1 | 0 0 0 || 0 0 0 | 0 0 -2 | 0 0 0 || 0 0 0 | 0 0 1 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 1 0 0 || 0 0 0 | 0 0 0 |-2 0 0 || 0 0 0 | 0 0 0 | 1 0 0 %
% 0 0 0 | 0 0 0 | 0 1 0 || 0 0 0 | 0 0 0 | 0 -2 0 || 0 0 0 | 0 0 0 | 0 1 0 %
% 0 0 0 | 0 0 0 | 0 0 1 || 0 0 0 | 0 0 0 | 0 0 -2 || 0 0 0 | 0 0 0 | 0 0 1 %
% =========================================================================================================== %
% 0 0 0 | 0 0 0 | 0 0 0 || 1 0 0 | 0 0 0 | 0 0 0 ||-2 0 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 1 0 | 0 0 0 | 0 0 0 || 0 -2 0 | 0 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 1 | 0 0 0 | 0 0 0 || 0 0 -2 | 0 0 0 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 1 0 0 | 0 0 0 || 0 0 0 |-2 0 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 1 0 | 0 0 0 || 0 0 0 | 0 -2 0 | 0 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 1 | 0 0 0 || 0 0 0 | 0 0 -2 | 0 0 0 %
% ----------------------------------------------------------------------------------------------------------- %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 1 0 0 || 0 0 0 | 0 0 0 |-2 0 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 1 0 || 0 0 0 | 0 0 0 | 0 -2 0 %
% 0 0 0 | 0 0 0 | 0 0 0 || 0 0 0 | 0 0 0 | 0 0 1 || 0 0 0 | 0 0 0 | 0 0 -2 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Axyz=ones(1,Nx*Ny*(Nz-1));
DZ2 = (-2)*diag(ones(1,Nx*Ny*Nz)) + (1)*diag(Axyz,-Nx*Ny) + (1)*diag(Axyz,Nx*Ny);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Building of the Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H=(-hbar^2/(2*m0*Mass)) * ( DX2/dx^2 + DY2/dy^2 + DZ2/dz^2 ) + diag(V0(:)*e) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% Diagonalisation of the Hamiltonien %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H=sparse(H);
[PSI,Energy] = eigs(H,n,'SM');
E = diag(Energy)/e;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% Normalization of the Wavefunction %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n
psi_temp=reshape(PSI(:,i),Ny,Nx,Nz);
%psi(:,:,:,i) = psi_temp / sqrt( trapz(z,trapz( y' , trapz(x,abs(psi_temp).^2 ,2) , 1 ),3) ); % normalisation of the wave function psi
psi(:,:,:,i) = psi_temp / max(psi_temp(:));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% here is a small patch due to differences between Octave and Matlab
% Matlab order the eigen values while Octave reverse it
if E(1)>E(2)
psi=psi(:,:,:,end:-1:1);
E=E(end:-1:1);
end
end